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ABSTRACT 


In  the  area  of  application  of  probabilistic  and  stochastic  modeling 
sequences  of  random  variables  evolving  over  time  are  usually  assumed  to 
be  sequences  of  independent  random  v^iables  or  Markov  sequences.  Here 
we  introduce  and  apply  a multivariate  Exponential''^  distribution  which  may 
describe  Markov  or  non-Markov  sequences.  The  present  work  has  examined 
one  particular  class  of  multivariate  exponential  distributions  which  preserve 
Markov  sequence  properties  for  both  modeling  of  downtime  distributions  and 
modeling  of  stages  of  component  deterioration.  In  downtime  modeling,  we 
study  the  distribution  of  the  sum  of  several  dependent  random  variables  and 
compare  the  result  with  the  distribution  of  a sum  of  independent  variables  as 
well  as  with  the  lognormal  distribution.  In  deterioration  modeling,  we 
consider  part  replacement  rules  based  on  observation  of  the  state  of  the  part's 
quality  and  on  specified  reward  structures.  We  identify  the  rate  of  deteriora- 
tion by  examining  how  long  the  component  stays  in  each  state  and  use  dynamic 
programming  to  set  up  recursive  optimization  equations  such  that  the  ex- 
pected reward  per  unit  time  is  maximized.  Sufficient  conditions  are  given 
under  which  the  optimum  replacement  rule  has  a very  simple  structure. 
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A n s T n ACT 

In  the  area  of  application  of  probabilistic  and  stochastic  modeling, 
sequences  of  random  variables  evolving  over  time  are  usually  assumed 
to  be  sequences  of  independent  random  variables  or  Markov  sequences, 
j Here  we  introduce  and  apply  a multivariate  "exponential"  distribution 

I which  may  describe  Markov  or  non-Markov  sequences.  The  present 

work  has  examined  one  particular  class  of  multivariate  exponential 
distributions  which  preserve  Markov  sequence  properties  for  both 
i modeling  of  downtime  distributions  and  modeling  of  stages  of  component 

deterioration.  In  downtime  modeling,  we  study  the  distribution  of  the 
sum  of  several  dependent  random  variables  and  compare  the  result  with 
the  distribution  of  a sum  of  independent  variables  as  well  as  with  the 
lognormal  distribution.  In  deterioration  modeling,  we  consider  part 
replacement  rules  based  on  observation  of  the  state  of  the  part’s  quality 
and  on  specified  reward  structures.  We  identify  the  rate  of  deteriora- 
tion by  examining  how  long  the  component  stays  in  each  state  and  use 
dynamic  programming  to  set  up  recursive  optimization  equations  such 
that  the  expected  reward  per  unit  time  is  maximized.  Sufficient  condi- 
tions are  given  under  which  the  optimum  replacement  rule  has  a very 
simple  structure. 
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Two  important  considerations  in  the  design  and  the  evaluation  of 
systems  are  their  capability  to  return  to  normal  operation  after  severe 
deterioration  and  that  after  occurrence  of  failures.  This  study  develops 
and  analyzes  stochastic  models  of  systems  whose  components  can  be 
replaced  without  occurrence  of  failure  and  can  be  repaired  when  the 
failure  does  occur.  The  analysis  focuses  on  a model  in  which  system 
states  are  distinct  at  any  time;  the  system  is  operating  (up)  at  one  of 
the  different  grades  of  performance  or  it  has  failed  (down)  and  is  at  one 
of  the  different  stages  of  repair. 

With  the  advent  of  modern  technology,  system  states  can  be  ob- 
served with  a non-destructive  testing  (e.  g.  , acoustic  wave,  infrared, 

laser,  x-ray etc.  ) without  inducing  performance  degradation  or 

interrupting  repair.  Also  the  data  is  easy  to  analyze  on  a special  pur- 
pose minicomputer,  a microcomputer,  or  even  a microprocessor. 

One  of  the  great  challenges  which  engineering  faces  is  the  develop- 
ment of  large  and  complex  systems  for  both  commercial  and  military 
applications.  Outstanding  examples  of  such  systems  are:  nuclear  power 
plants,  sw  itching  machines  of  telephone  office,  air  traffic  control  sys- 
tems, radar  detection  and  warning  systems,  space  vehicles,  aircraft, 
communication  networks,  real-time  computer,  and  medical  instruments. 
Failures  of  these  complex  systems  might  be  fatal  to  human  life.  Re- 
placement of  deteriorating  parts  will  reduce  the  probability  of  failure. 
Once,  a failure  has  occurred,  repairability  is  an  important  considera- 
tion for  restoring  the  system  back  to  normal  as  soon  as  possible. 
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Tfif  111  ni,i jnti'nanci'  poli'  n's  for  niaintainabJc  and  repair- 

able systems  makes  use  of  information  about  prolnibility  distributions  of 
eomponent  lifetimes  and  downtimes.  Such  desij^ns  are  (aeilitalc'd  it 
these  distributions  have  simple  analytic  forms. 

Downtime  distributions  have  been  frequently  modeled  as  lognormal, 
Weibull.or  Erlang,  in  form  f 7 ].  These  distributions  are  all  skewed  and 
corresponding  to  non-negative  random  variables  like  downtimes.  Here 
we  consider  one  more  family  of  distributions  which  has  some  physical 
motivation. 

Since  a downtime  interval  is  often  the  sum  of  subsidiary  intervals 
(for  failure  isolation,  component  reinoval,  repair,  reassembly,  align- 
ment, etc.  ) it  seems  reasonable  to  think  of  the  downtime  as  the  sum 
of  subsidiary  time  intervals: 


n 


The  subscript  of  x reminds  us  of  the  number  of  summands.  Several 
' n 

distributions  are  possible  for  the  individual  r but  here  we  consider 
exponential  distributions  which  are  the  simplest  and  are  also  widely  used 
to  represent  random  times  between  events.  Clearly,  once  the  distribu- 
tion of  such  an  x has  been  characterized,  the  assumption  that  r is 
n ' 

exponentially  distributed  can  l)e  weakened  to  allow  an  x^-typ<'  distribu- 
tion for  any  r^.  Eor  example,  the  repair  time  interval  may  not  be 
characterized  by  an  exponential  distribution.  However,  by  using  the 
above  argument,  we  can  find  a suitable  n for  an  x^-type  distribution 
for  the  repair  time. 
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It  is  v,fll  known  that  the  downtime  will  have  an  I'.rlang  distri- 
bution il  the  r «ire  i nclep*- ncle nt  t'xjonential  random  variables  with 
1 

identical  mean  values.  Muth  [22]  has  considered  the  approxim<.ition  of 
Weibull  and  lognormal  distributions  by  x^  in  which  the  r^  are  indepen- 
dent exponential  variables  but  with  possibly  different  mean  values.  Here 
we  generalize  his  w'ork  to  allow  dependence  among  the  a rea- 

sonable situation  if  the  variables  represent  related  steps  in  a sequence 
of  downtime  operations  or  in  a sequence  of  deterioration  levels. 

Several  types  of  multivariate  exponential  distributions  have  been 
proposed  for  various  reliability  applications  [ 2 ],  The  multivariate 
exponential  distribution  discussed  most  is  the  Marshall  and  Olkin's  shock 
model,  Marshall  and  Olkin's  motivation  was  mainly  for  a system  with 
parallel  time  intervals  which  happen  to  start  at  the  same  time  (Figure 
1(a)). 


■*  r 
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(a)  (b) 

Figure  1: 

Parallel  and  Series  Random  Time  Intervale  . 


Time 


For  our  purpose  we  would  like  to  have  a sequence  of  random  tim«‘  inter- 
vals that  happen  one  by  one  as  in  Figure  1(b). 

The  present  work  has  examined  one  particular  class  of  multivar- 
iate exponential  distribution  for  both  modeling  of  downtimes  and  modeling 
of  component  deterioration. 


r 
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Chapter  II  (1 1 sfus  SI- s srvoral  I'ornis  of  imilti  var  i..U'  <1  i si  r i htil  ions 
which  have  exponential  marginal  distributions.  The  one  based  on  the  sums 
of  the  squares  of  normal  variables  is  selected  because  of  its  analytical 
s implic  ity. 

Chapter  III  gives  examples  of  the  distributions  of  the  sums  of  n 
dependent  exponential  variables.  We  show  that  the  introduction  of  de- 
pendence among  r^^  (with  possibly  unequal  means)  does  not  broaden  the 

class  of  X distributions  over  that  which  results  from  independent  r.. 

n 1 

That  is,  the  sum  of  n dependent  exponential  variables  has  a distribution 
identical  to  that  of  the  sum  of  n other  independent  exponential  variables. 
Several  conclusions  are  on  the  approximation  of  lognormal  variables  by 
sums  of  exponential  ones,  along  with  other  possible  extensions. 

Chapter  IV  examines  maintenance  policies  for  systems  in  which 
the  degree  of  deterioration  can  be  observed  continuously,  A Markov 
sequence  model  is  developed  where  the  holding  times  in  the  varmus 
states  are  multivariate  exponentially  distributed,  V'e  assume  that  the 
functioning  rewards  during  the  system's  lifetime  decrease  with  the  in- 
creasing deterioration.  Additional  dependence  properties  of  multivariate 
exponential  variables  are  developed  for  the  optimization  studies.  The 
dependence  relations  seem  to  be  consistent  with  what  one  expects  of  a 
general  system  in  the  real  world.  Optimal  replacement  rules  which  maximize 
the  expected  reward  per  unit  time  are  obtained  by  using  the  dynamic  program- 
ming method.  Sufficient  conditions  which  simplify  the  optimal  decision 
rules  are  given  in  Appendix  1.  Chapter  V summarizes  the  results  and 
makes  suggestions  for  further  research. 

Part  of  Chapter  II  and  111  have  already  appeared  in  a previous  report 
by  Hsu  and  Shaw  [ IZ]  . 

.. 
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MUI/J  IVAKIATK  I- XPONKNTIAI.  DLSTlUBU'i  IONS 

The  mul tinormal  distribution  has  been  studied  far  more  exten- 
sively than  any  other  multivariate  distribution.  Indeed,  its  position  of 
pre-eminence  among  multivariate  continuous  distributions  is  more  marked 
than  that  of  the  normal  among  univariate  continuous  distributions.  How- 
ever, in  practical  engineering  problems  there  have  been  signs  that  the 
need  for  multivariate  exponential  distributions  is  becoming  recognized. 
Several  forms  of  multivariate  exponential  distributions  have  been  intro- 
duced, The  one  based  on  the  sums  of  squares  of  multinormal  variables 
is  selected  here,  because  of  its  analytical  simplicity  through  transfor- 
mation to  well  developed  multinormal  distributions. 

This  chapter  reviews  several  popular  multivariate  exponential 
distributions  and  develops  properties  of  the  one  to  be  used  in  later  chapters. 
Some  important  dependence  properties  will  be  developed  in  Chapter  IV. 

2,  1 Review  of  Literature 

It  is  well  known  that  a multivariate  distribution  is  not  uniquely 
specified  by  its  marginal  distributions  [ 2 ].  For  example,  while  the 
bivariate  normal  distribution  has  nice  analytical  properties,  it  is  not  the 
only  one  with  normal  marginals.  This  multipli'-ity  will  be  demonstrated 
for  the  case  of  exponential  marginals  by  considering  a few  possible  bi- 
variate densities. 

Gumbel  [ 9 ] considered  several  bivariate  exponential  densities. 

The  first  Fj(r^,  r^)  is  based  on  the  following  general  formula  for  com- 
bining marginal  distributions; 


-D- 


F(x.  y)  = Fy(x)  FY(y)  fl  + a[  1-F^(x)]r  l-Fyfy)]  } . ( 2.  1 ) 


O'  I < 1 


The  marginal  distributions  of  X and  Y are  each  exponential.  When  ap- 
plied to  exponential  variables  with  unit  mean  values  this  produces  the 


distribution  function 


F^(r^,  r^)  = (1-e  ^)(1+Q'e 


(2.  2) 


rj.  r^  > 0 . 


O'  < 1 , 


and  the  density  function: 


r^)=  e 


[1  +a(2e  *-l)(2e  ^-1  )] 


(2.  3) 


In  this  model,  ct  = 0 corresponds  to  independence  of  r^  and  r^,  and  it 


can  be  shown  that  the  correlation  coefficient  p 

^1*  *^2 


*'2 


= 0/4 


(2.  4) 


with  its  magnitude  limited  to  be  less  than  or  equal  to  l/4.  Another  model 
of  Gumbel*s  is  defined  by  the  joint  distribution  function: 


"•l  -■■2  -'■|-'■2-« 

f-e  -e 


0<G£l  , — ® 


and  the  corresponding  density: 

- r. - Ty-Or  r^ 

f2(r^.r2)=e  ’ ' ^ [ ( t + 0 r ^ )(1  + ©r  1 


(2.6) 
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) It  I I f/  ii  i ( 1 1 I I |i'  * * I ^ *>  * *^**  ^*'  P**  ^*^^  * ’ ^*^  * * a/it'l  . 

p = E (e"S-l  . 

‘•r'-z 

where:  ^ 

E.(6'''  ) = I Z"^dZ  . 

0-^ 

The  univariate  exponential  distribution  derives  considerable  im- 
portance from  its  role  as  the  distribution  of  waiting  time  in  a Poisson 
process.  It  is  natural  to  inquire  whether  a similar  relationship  exists 
between  some  bivariate  exponential  distribution  and  the  waiting  times  in 
a suitably  defined  tv^’o-dimensional  Poisson  process.  One  such  possi- 
bility, investigated  by  Marshall  and  Olkin  [ 2 ] and  later  generalized  by 

them  [ 21  1 . will  now  be  described. 

This  distribution  can  be  thought  of  as  a result  of  fatal  shocks 
occurring  from  three  independent  Poisson  sources  with  rates 

Component  1 with  lifetime  r^  is  killed  by  events  of  the  first  or 
third  variety,  and  r^  is  determined  by  events  of  the  second  or  third 
type.  We  can  define  the  reliability,  which  is  a probability  of  survival  as: 


R{r,,  rp  = P [ R,  ^ r,,  R^  > r^l 


or 


F3(r,.r^)=  1 - 


.(A,+X|3)r| 


+ e 


-riXj-r2X2-X^2max(rj,  r^) 


Here  the  correlation  coefficient  is  : 


(2.7) 


(2.  8) 
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C'l.  ‘n 


which  ranges  bet'ween  zero  and  one. 

Barlow  and  Proschan[  2]  point  out  that  Eq.  (2.8)  is  the  unique  bi- 
variate exponcntialdistribution  with  the  zero-memory  property;  the 
joint  survival  probability  of  a pair  of  components  each  of  age  t is  the 
same  for  all  t (e.  g.  , the  same  as  if  both  were  new).  This  zero  memory 
is  quite  desirable  when  modeling  joint  lifetimes  of  components  in  a sys- 
tem. In  the  present  context  of  r.  representing  durations  of  a sequence 
of  related  events,  this  special  property  seems  inessential. 

Kibble  [ 1 b]  considered  a bivariate  exponential  density  of  the 

form  : 


’■l^"2 

/ 2 

1 

“ ^ 2,,  2. 
2o-  (1-p  ) 

^0 

4<r'’(l-p^) 

cr^(l-p^) 

^ « 

- 

in  w 


hich  I (•  ) is  a modified  Bessel  function  of  order  zero  defined  as  ; 


(z/2)^^ 


k=0  (k-L) 


(2.  11) 


A detailed  study  of  this  distribution  has  been  given  by  Nagao  and 
Kadoya  | 23  1 • 

2.2  Multivariate  Exponential  Distribution 

We  will  derive  Equation  (2.10)  in  an  alternate,  more  natural  way 
for  use  below.  The  same  approach  has  also  been  taken  by  other  authors, 
for  example,  Moran  and  Vere  - Jones  [ 32]  . 
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Thf  fU  iihily  ill  I j!.  10)  apyjlii-s  v«/h<-ii  w*-  vi<'v<  anil  as  (x-inn 

y'ene  rated  froni  correlated  normal  variables.  It  is  well  known  that  if  w^ 

2 

and  Zj  are  independent,  zero  mean,  equal  variance  (it  ) normal  vari- 
bles  then  r^  defined  as  ; 


= w.  + z, 


(2.  1 2) 


hcas  an  exponential  distribution  with  mean  : 

E(r  ) = . (2.  13) 

r 

Now  if  (w^,  w^)  and  (z^,  z^)  are  two  independent  pairs  of  normal  vari- 
able s,  but  with : 


covfw.,  w.)  = cov(z.,  z.) 
1 J 1 J 


i = 1,  2:  j = 1,  2 


( 2.  1 4) 


then  r^  and  r^  defined  by  : 


2 2 • 4 , 
r.  = w.+z.  1=1,2 

111 


( 2.  1 5) 


will  be  dependent  exponential  variables. 

Equation  (2.10)  can  be  derived  from  this  reasoning  when  all  four 

2 

variables  (w^,  w^,  z^,  z^)  have  equal  variances  ir  . In  that  case  : 


f(Wj.W2,z^,Z2)  = 


(2  7r)^(1-p^)(r^ 


exp 


Wj  -2pWj w^-hw^  + z^-2pZjZ2+Z2 


2(l-p^)(r^ 


and  : 


( 2.  1 h) 
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2^2., 
Wi  + r, 


2^2,^ 
W2  + ^2  - 2 


( 2.  1 7) 


Introduction  of  polar  coordinates  in  w^,  and  w^.  planes: 

cos  . "^2  ~ '^2.  ^2 

Zj  =7,  sin  0 ^ = 72  ^2  ’ 


reduces  the  integral  to  the  form 


s/r,  \'r 


1 1 ''.'2  27r  2;r 

^2^  = “TT — f ' 

47r  (1-p  )o-  ^^  = 0 7^=0  ej  = 0 02=° 


7j  +72  “ 2071^2 


)T 


• 7,72  S7|d72d0^d02 


( 2.  1 H) 


Substitution  of  cp  = 0,-02  °1  integration  produces  aperiodic 

integrand,  independent  of  02-  Thus  the  0 -integrals  become: 


2 n- 

Zn 

Z-n 

f 

f - 

■■  Zn  J 

o 

tl 

02=0 

0 

.2/,  „2, 


dtp 


47r^  F 


Ip  I 7,^2 


(2.  19) 


-It  - 


/■■ina}ly,.s.ib«lihj1ion  of  (2.  19)  into  (2.  IK)  and  differentiations  with 
respect  to  r^  and  r^  produces  (2.10).  It  turns  out  that  is  just 

the  square  (p^)  of  the  correlation  of  the  underlying  normal  variables,  with 

a range  between  zero  and  one. 

The  Bessel  function  form  for  this  bivariate  exponential  based  on 
normal  variates  may  not  appear  to  be  very  felicitous.  However,  this 
kind  of  distribution  will  be  convenient  when  we  concentrate  in  the  next 
chapter  on  the  sum  of  dependent  exponential  variables,  as  in  ( 1.  1 ). 

For  the  n-dimensional  version  of  this  class  of  distributions,  we 
consider  zero  mean,  normal  n-vectors  w and  z each  with  the  same  co- 
variance  matrix  (which  is  positively  definite  and  symmetric)  : 


E[  w w']  = E[^_z']  = r . 


( 2.  20) 


In  this  way,  for  each  i,  w.  and  z.  will  have  the  same  variances  so  the 
sum  of  their  squares  will  be  an  exponential  random  variable  r..  We  do 
allow  r.  and  r^  to  have  unequal  means,  contrary  to  the  special  case 
in  (2.  10).  Using  the  underlying  normal  distributions  it  is  easy  to  show' 


that  the  r^  have  means  ; 


El  r,l  = , 


(2.  21) 


and  correlation  coefficients 


(V  /‘■'iiV  i ■ 

1’  J 


(2.  22) 


which  can  take  on  values  between  0 and  1. 

The  general  approach  to  calculation  of  the  multivariate  n- 

distribution  is  through  integration  of  the  normal  density  Kw^.w^, 

^2'  • • • * appropriate  regions,. 


A mullivui-ial.-  < abe  of  interest  corresponds  to  the  covariance 


iii.ii  r IX  : 


•v  = (T  T.p 
II  I .1 


l-I 


( 2.  22) 


in  ^huh  correlation  between  r.  and  r.  falls  off  exponentially  in  relation 
to  the  separation  between  their  indices  i and  j.  The  corresponding  joint 
density  has  a nice  structure,  revealed  by  its  trivariate  form; 


1 

f(r  . r^.  r ) = ? 2 ZT.  ^ 


is 


I 


i 


^ 0-  0-  (1  -p^) 

I 1 ^ 


• ^0 

/ 2 

/.  2. 
O-^O-gd-p  ) 

J 

-1 

(2.  24 ) 


This  expression  generalizes  in  an  obvious  way  to  the  n-variate  form: 

■ 1 n-1 


/ I 2 

N/r.  r. 


n-l  n -1  n-t 

= 


(i-p^), 


1 


^ 2 *'n 


exp  - 

1 2(t-p‘')\(r^  i=2  or": 


(n  ^ 2).  (2.  25) 

The  special  case  of  Eq.  (2.23)  describes  a Markovian  normal  sequence 
w w ...  w with  stationary  correlation  coefficients.  The  sequence 

I’  2’  n 

r j,  r^,  . . . r with  a joint  density  of  (2.25)  can  be  easily  seen  to  be 
Markovian  also,  a fact  also  established  in  Griffiths  [ IB]  . 
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CH  APT  Eli  _ I IJ 
IJOWN'IIMJ',  MOlJi  lANCj 

As  mentioned  in  Chapter  I,  we  want  to  represent  a downtime  as 
the  sum  of  several,  possibly  dependent,  subsidiary  times  : 
n 

X = S r . (3.  I ) 

^ i=l  ^ 

We  will  examine  several  possibilities  for  r^^  distributions,  beginning 
with  a review  of  Muth' s work  [ 22]  for  the  case  of  independent,  expo- 
nentially distributed  r^  with  unequal  mean  values.  Then  we  will  gene- 
ralize his  model  to  allow  dependence  among  the  r^^,  seeing  that  the  family 

of  possible  f is  not  expanded  in  this  way,  but  that  it  allows  for  con- 
n 

venient  study  of  the  case  of  large  n and  possible  convergence  to  a log- 
normal density. 

3.1  Review  of  Muth*  s Results 

In  Appendix  B of  Muth' s thesis  [ 22],  he  investigated  the  proper- 
ties of  the  class  of  probability  distributions  belonging  to  the  family  of 
Erlang  distributions,  which  are  generated  as  the  distribution  of  a sum 
of  independent  and  exponentially  distributed  random  variables.  Speci- 
fically, he  compared  these  distributions  with  the  gamma,  lognormal, 
and  Weibull  distributions,  which  are  well  known  examples  of  two- 
parameter  unimodal  distributions  and  which  have  been  used  as  models 
for  repair  times. 

He  considered  the  random  variable  Y(n)  which  is  the  sum  of  n 
independent  random  variables  X.  : 

Y(n)  = f . . . + , (3,  2) 


-1-U 


1 


whert!  (!ach  X.  is  e^xpunf  nl  ial  ly  (list  ri  lulled  wifli  nu-.iii  value  iiu.  I lu 
raiidoui  v.iriahli  Y(ii)  ii.iH  a d i s ( ri  lui  I i oti  in  wliiili  n p.t  r.,  i nt- 1 e rn,  ilu 

values  of  m to  m , must  be  specified,  and  where  n itseli  is  a vai  iable. 
In 

We  use  (he  following  notations  and  definitions; 


Et  Y(n)]  = m. 


Var(  Y(n)l  = 


2 

0-  , 


(3.  3) 


skewness)  = 


y (skcwil^oij/  — Q 


(3.  4) 


EUXtoU  . 3 . 
0* 


(3.  5) 


In  order  to  compare  the  probability  law  of  Y(n)  with  that  of  a random 
variable  Z,  having  a two  parameter  distribution,  he  set  their  means 
and  variances  equal : 


E[  Y(n)]  = E[Z1  . 

(3.6) 

Var  [ Y (n)]  = Var  [ Z] , 

(3.7) 

and  compared  the  Iiigher  moments  of  both  probability  laws.  He  re- 
stricted the  comparison  to  the  third  and  fourth  moments.  The  under- 
lying idea  is  that  the  approximation  of  one  probability  law  by  another 
becomes  closer  as  one  successively  matches  higher  moments,  and  that 
matcliing  the  third  moment,  together  with  the  constraints  (3.6  ) and 
( 3,7  ),  does  already  provide  an  improved  approximation.  Without  loss 
of  generality  all  distributions  under  consideration  were  normalized, 
such  that  their  mean  values  equaled  1.  The  relative  dispersion  of  a 
distribution  is  expressed  by  its  coefficient  of  variation  y , namely  for 
any  Z : 


I 

I 

( 


1 


i 


-15- 


) 


^/Va  r(Z) 

.M/) 


Cl.  H) 


He  studied  the  ranges  of  7^,  7^  and  v with  the  above  constraints 
on  the  mean  and  variance  of  Y(n)  and  found  that  those  parameters  were 
maximized  by  finding  two  numbers  m^  and  mj^  and  setting: 


m,  = m , 
1 a 


m.  = m,  , 
1 b 


i = 2,  3. 


(3.  9) 


Similarly,  those  parameters  were  minimized  by  setting: 


m.  = m i = 1,  2,  . . . , n-1, 

X cL 


m = m,  , 
n b 


(3.  10) 


for  the  same  m^  and  If  v is  fixed,  then  7^  and  7^  are 

according  to 

2 ■>  7|  (”)  ^ " I ^ > 0, 

6 '>  y 2}''^'^  >2  — ^ 


bounded 


(3.  11) 


Increasing  n further  increases  these  maximum  values  7^  and 
whose  limits,  for  n -*  are:  max  7^  = 2 and  max  = 6.  Values  of 
^ and  72^®  functions  of  n were  computed  for  the  gamma,  lognormal, 
and  Weibull  distributions,  and  are  presented  in  Figures  2 and  3 together 
with  the  feasible  Muth' s generalized  Erlang  distributions.  These  graphs 
show  that  the  Erlang  distribution  can  be  used  to  approach  the  lognormal 
distribution  in  the  third  and  fourth  moment,  if  the  coefficient  of  variation 
is  less  than  0.  6.  The  Weibull  distribution  on  the  other  hand  cannot  be 
approximated  in  this  fashion. 


0 0.2  0.4  0.6  0.8  1.0 


Coefficient  of  Variation  v -* 

Figure  2:  Coefficients  of  Skewness  (from  E.  J,  Muth, 
Repairable  Systems,  Ph.  D.  Dissertation, 
PIB,  1967)  . 


Figure  ; Coefficients  of  Kurtosis  (from  E.  J.  Muth, 
Repairable  Systems,  Ph.  D.  Dissertation, 
PIB.  1967)  . 
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n.  2 Sums  of  Correlated  Exponential  Variables 

'J  hf  characteristic  function  of  x defined  in  ( 1.1  ) as  the  sum  of 

n 

exponential  random  variables  can  be  computed  iis  follows,  Usin^  (2.  1 h) 
and  the  definitions  : 


"2  *^2 

V = £ w y = £ z.  , 

1 1 


(3.  1 2) 


je  can  write  the  characteristic  functions  as  : 


(s)  = 4>  (s)  fl*  (s)  = ^y(s)  , 

^n  V y 


due  to  the  independence  and  identical  distributions  of  w and  _z.  The  possi- 
bly correlated  variables  w^  can  be  represented  as  linear  transformations 
of  independent  unit  variance  normal  variables 


w = MiC,  Ef^£'  ] = I 


(3.  13) 


(Here,  as  throughout  this  report  I is  the  identity  matrix  of  the  appropriate 
size.  ) Similarly  : 


z = Mi,  Elii'l  = I 


(3.  I'll 


It  follows  that 


V = w'  w = 1C'  M'  M C . 


(3.  15) 


<1)  (s)  = (^/2^)  (...  r exp^-  - si'M'Mcl  rf? 

V .}  J ^ j 

= (vT^)  f. . . / exp  i i'  R“^il  di  , 


(3.  lb) 
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whrrt;  we  liave  defined  Ihe  mafrix  : 


- 1 


1 


lf>sM'M 


(3.  17) 


The  integral  in  (3.  1 6)  is  of  the  form  of  a normal  density  integrated  over 
all  values,  except  for  a scale  factor.  Thus: 

-1 


.l>^(s) 


= (JJ^ 


(3.  1H) 


and  the  desired  characteristic  function  for  the  sum  of  exponential  vari- 
ables is  : 


<))  (s)  = ( )l  + 2s  M'  M|  ) = l/q(s)  . 


-1 


(3.  19) 


Equation  (3,  19)  shows  that  the  characteristic  function  of  is 
the  reciprocal  of  an  n -degree  (or  less,  the  rank  of  T is  less  than  or 
equal  to  n)  polynomial  q(s)  whose  coefficients  are  determined  by  the 
covariance  matrix  F.  Properties  of  that  polynomial  characterize  x^. 

In  particular,  we  have  the  theorem; 

n 

All  possible  density  functions  f for  x^  = T r^  can  be  achieved 

n 1 

with  independent  r^,  i.  e,  , with  diagonal  F. 

The  proof  begins  by  noting  that  the  roots  of  q(s)  are  negative  re- 
ciprocals of  the  eigenvalues  of  the  symmetric  nonnegative  definite, 
matrix  ZM'M  (or  equivalently  of  2 F),  so  they  are  negative  numbers.  If 
the  r,  are  independent,  then  M is  diagonal: 


Mj  = diaglo  j . . . , a^) 


(3.  20) 


and  so  is  : 
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/ 2 2 

2M^  = flia)4(2o'  ^ • • • . Z'-"  ■ (:■!.  21  ) 

In  this  "inch-pendent"  ease  : 

q (s)  = ^ (1  + .s)  • (3.  22) 

^ i = l 

The  prooi  is  completed  b-y  comparing  the  pol-ynomial  q(s)  in 
(3.  19)  for  a general  M with  the  qj(s)  in  (3.  22)  for  independent  r..  The 
mean  values  of  the  latter  (2a  can  be  chosen  to  make  the  roots  of  qj.{s) 
match  an-y  n (necessarily  real)  roots  of  q(s).  The  resulting  polynomials 
and  characteristic  functions  will  be  identical  because  the  necessary 
4>{0)  - 1 property  removes  any  scale  factor  ambiguity.  This  completes 
the  proof.  I1 

In  general,  the  n-independent  exponential  variables  whose  sum 
is  indistinguishable  from  the  sum  of  n-correlated  ones  will  have  differ- 
ent mean  values  from  those  of  the  correlated  variables. 

1.  3 Large  Number  of  Summands  and  Lognormal  Approximations 

Once  this  structure  has  been  established  for  sums  of  correlated 
exponential  variables,  previous  results  for  stuns  of  independent  variables 
(e.  g.  , those  of  Muth  [ 22)  ) are  directly  applicable.  He  sought,  among 
other  results,  appropriate  mean  values  for  independent  exponential  var- 
iables such  that  their  sum  should  well  approximate  a lognormal  variable 
whose  density  is  of  the  form: 

f(x)  = (x.r  n/^)  t 2<t'  , Q _ (3.23) 

with  mean  and  variance: 
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I 

I 


r?  f o-  “/2 


V a r f x] 


2f)  f I) 


:m) 

(;{.  21-)) 


While  the  theorem  above  suggests  there  might  be  no  advantage 
to  allowing  correlations  among  the  r^^  when  the  in  ( 1 , 1 ) is  to  ap- 
proximate X,  experience  has  shown  that  this  new  viewpoint  can  be  con- 
venient. Figure  4 shows  a lognormal  density  having  Fix]  = 1 and 
Varfx]  = 0,  653646,  along  with  approximations  to  it. 

The  approximations  in  Figure  4 are  sums  of  correlated  expo- 
nential variables  based  on  underlying  normal  distributions  having  covar- 
iances of  the  form  of  (2.  23).  Fur  the  rmore.the  sunimands  are  assumed 
to  have  equal  mean  values  (all  cr^  = cr| ).  In  this  way  each  x^  is  com- 
pletely characterized  by  three  numbers;  n,  p,  (t^. 

The  mean  values  are  chosen  to  be  : 


Ffr.]  = 2(r^  = ^ , (3.  26) 

so  the  mean  of  x matches  that  of  x.  The  correlation  parameter  n in 
n ‘ 

( 2.  23)  is  then  chosen  so  that  the  variances  of  and  x are  equal.  It  is 
straightforward  to  compute: 


var(x^)  = (1  -h  p ^)/2  , 

var{x^)  = (3  + 4p“  + 2p^)/9  , 

var(x^)  = (2  + 3p^  + 2p^  + p^) /8  . . (3.27) 


using  the  underlying  normal  densities  and  the-  mean  values  de-terniined  in 
( i.  26).  For  the  specified  Var(x)  = 0.  653646  it  turns  out  (hat  the-  e orre-- 
lation  factors  must  be  as  shown  in  Table  t. 
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I.I  1 

( .1  > r 1' 1- 1-it  ii 111  I'  .ii  liir»  lor  No r 1 1 i.i  1 1 /.ei  1 

n 

n 

p 

V 

0.  554339 

3 

0.  750000 

! ^ 

0.  820321 

In  principle,  the  densities  can  be  computed  by  Laplace  trans 

n 

form  inversion  of  { 3.  1 9)  after  the  p and  (t.  values  have  b«.en  specified. 
An  alternative  approach  is  possible  once  the  roots  (-X.)  of  q(s)  have 
been  determined.  For  example,  if  all  those  roots  ar<;  distinct  then: 


n 


n 1=  1 


a. 

1 


e 


-X.t 

X 


C^.  28) 


One  condition  on  the  unknown  is  : 

f f (t)dt  = 1 = S a./X,.  . 

' X 11 

on  1 

Application  of  the  Initial  Value  Theorem  to  (3.  28)  shows  that  derivatives 
of  the  density  must  satisfy  : 


f^^^  (0)  = 0 

X 

n 

n , 

= 2 (-X.)  a.;  k = 0,  1,  . . . , (n-2)  . 

^ i 1 


(3.  30) 


P.'quations  (3.  29)  and  (3.  30)  provide  n-equations  for  the  n-unknowns  a.  . 

This  approach  has  been  used  to  compute  the  approximating  den- 
sities, with  the  following  results:  (single  precision) 
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, = ,.»n3<.r.  ,,-l.  2X672.  48772, 

X 


1 = 0.  99ir)90  + 1.67508  e 

^3 


- 2.  6b67  e 


-6.  8571  41 


I - 1.65610  e'*-  26572.  _ (,9305  5-2-79522. 


6646  .-2'- -0.  145528  e-X-'-'’-'*"*  . «•  31 ) 


Figure  4 reveals  that  as  n increases,  this  particular  kind  of 
sum  of  correlated  exponential  variables  has  a distribution  which  ap- 
proaches the  lognormal  shape.  Table  2 shows  the  mean  values  for  the 
equivalent  independent  exponential  random  variables  whose  sums  have 
the  same  distributions,  respectively.  These  means  are  the  reciprocals 
of  the  exponents  in  (3.  31).  There  is  considerable  intuitive  appeal  to  this 
kind  of  a limit  of  summing  more  and  more  terms,  each  contributing 
equally,  on  the  average,  with  correlation  between  two  summands  de- 
creasing exponentially  as  their  separation  increases.  However,  no  pre- 
cise asymptotic  result  seems  possible.  The  " heavy- tail"  of  the  log- 
normal (slower  than  exponential  decay  of  f(x)  for  large  x)  wil  never  be 

matched  by  an  x since,  for  every  n,  the  tail  of  f^  (x)  will  be  domi- 

n 

nated  by  the  slowest  decay  of  its  (at  most  n)  decaying  exponential  sum- 
mands. 

Other  obstacles  to  having  the  distribution  of  a sum  of  exponential 
variables  approach  the  lognormal  form  are  given  in  Muth's  analysis  of 
bounds  on  the  moments  of  such  sums.  Those  results  for  independent 
summands  are  applicable  to  correlated  ones,  due  to  the  theorem  pre- 


sented earlier. 
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lAlil-l  2 

Means  of  Independent  Exponential  Variables  Whose  Sums  Have 

Densities  f 

X 


n 

— 

^2 

Bi 

’^4 

2 

0.77716 

0.  2228 

3 

0.  79286 

0.  1458 

0.  061 3 

i ^ 

0. 79636 

0.  1 2828 



0.  04673 

0.  02863 

The  foregoing  example  was  based  on  a covariance  matrix  of  the 
form  of  (?.  23)  in  which  o and  p were  chosen  so  x^  would  have  pre- 
scribed mean  and  variance.  This  approach  could  be  generalized  to  co- 
variance  matrices  having  more  adjustable  parameters  in  the  hopes  of 
matching  more  moments  of  either  experimental  data  or  of  the  lognormal 
distribution.  One  example  would  be  the  following  with  parameters  cr , 


7ij= 


li-J 


+-  Bt:  , 


D-J 


) 


(3.  32) 


whe  re 


(TT,  4 7r^)‘ 


— ^ - |7rj+7r,^| 


“ ^2^ 

(tT  J • 
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lliivwi'i  I,  Millli'rt  1)01111(1  on  lln  h 1' *■  w lie  K i’  ol  .1  Iiiilti  ol  (•  x J)ol)<' ri(  i oi  s 
precludes  adjustment  ol  the  three  parameters  here  to  mateli  (he  lirst 
tliree  moments  ol  every  possible  distribution. 

3.  4 Conclusions 

A multivariate  exponential  distribution  based  on  sums  of  squares 
of  normal  variables  was  examined.  We  showed  that  variation  in  the  cor- 
relation and  means  of  such  variables  does  not  produce  any  distributions 
for  their  sum  which  could  not  have  been  realized  by  a sum  of  indept'ndent 
exponential  variables.  However,  sums  of  such  variables,  when  they  have 
equal  means  and  stationary  correlations,  tend  to  have  distributions 
shaped  like  the  lognormal,  except  in  their  tail  regions. 

It  should  be  noted  that  much  of  the  popularity  that  the  lognormal 
distribution  enjoys  for  downtime  modeling  probably  results  from  the 
easy  use  of  normal  probability  paper  for  estimation  of  lognormal  para- 
meters from  data.  In  most  cases,  other  distributions,  say  f^  for  some 

n 

small  n and  suitable  <t ^ and  p,  might  fit  the  data  equally  well.  In  any 
«‘vent,  the  peculiar  lognormal  tail  behavior,  corresponding  to  rare 
events,  is  not  likely  to  be  a necessary  property  of  a distribution  being 
fitted  to  typical  amounts  of  data. 
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X 

Figure  -i;  Lognormal  Density  and  Approximations  to  it 

mean  = 1 

variance  = 0,  653646 

coefficient  of  variations  0.  808484 
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C H A P T P K IV 

OPTIMAL  HI  PLACKMKNT  OF  DI’TFRIOHA'J  ENG  PARIS 

A siiiiili-  (. omjxjin  nt  system  is  assumed  lu  protiiess  (tirou^h  a 
finite  number  ni  increasingly  bad  levels  of  deterioration.  The  system 
with  level  i(0  f_  i 5.  n)  starts  in  state  0 when  new,  and  is  definitely 
I replaced  upon  reaching  the  worthless  state  n . It  is  assumed  that 

the  transition  times  are  directly  monitored  and  the  admissible  class  of 
strategies  allows  substitution  of  a new  component  only  at  such  transition 
times.  The  durations  in  various  deterioration  levels  are  dependent  ran- 
dom variables  with  exponential  marginal  distributions  and  a particularly 
convenient  joint  distribution.  Strategies  are  claosen  to  maximize  the  aver- 
; age  rewards  per  unit  time.  For  some  reward  functions  (with  the  re- 
ward rate  depending  on  the  state  and  the  duration  in  this  state)  the  i 

knowledge  of  previous  state  duration  provides  useful  information  about 
the  rate  of  deterioration. 

I ‘ 

4.  1 Review  of  Literature  and  Introductory  Remarks 
i Many  mathematical  studies  have  been  devoted  to  the  optimization 

of  rules  for  inspecting  system  quality  and  for  repairing  or  replacing 
j parts  as  they  are  observed  to  deteriora'e.  This  section  reviews  some 

j of  those  previous  results  and  proposes  a new  model  which  allows  use 

' of  a measure  of  deterioration  rate  by  the  controller  which  replaces 

parts  so  as  to  optimize  the  av'erage  reward  per  unit  time. 

1 

' Many  inspection  policy  models,  where  inspections  reveal  a mal- 

I 

; function  occurrence  only,  have  been  published;  for  example,  Luss  and 

Kander  fl9].  Several  papers  deal  with  models  where  deterioration  can 
' be  observed,  and  most  of  them  concentrate  on  replacement  policies 

i 

1 

i 
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wliic  ti  assuini'  lliat  the  system's  state  is  always  known.  Herman  f -I], 

( r,  |,  Harlow  and  Proschan  f 1 |,  and  others  studied  such  models  as- 
suming that  the  deterioration  process  is  described  by  a transition 
probability  matrix  of  a Markov  chain.  Kao  [i6]  studied  optimal  re- 
placement rules  when  changes  of  states  are  semi-Markovian.  Rosen- 
field  [26]  examined  properties  of  optimal  policies  for  models  in  which 
the  system's  state  is  observed  by  inspections  only.  Kander  [15]  also 
examined  inspection  models.  However,  he  assumed  that  the  operating 
costs  occurring  during  the  system's  life  do  not  change  with  the  increas- 
ing deterioration. 

Luss  [20]  also  examined  inspection  models,  he  assumed  that  the 
operating  costs  occurring  during  the  system's  life  increase  with  the 
increasing  deterioration.  However,  he  assumed  that  the  holding  times 
in  the  various  states  are  independently,  identically,  and  eicponentially 
distributed.  The  policies  examined  include  the  scheduling  of  inspec- 
tions (when  an  inspection  reveals  that  the  state  of  the  system  is  better 
than  certain  critical  state  k)  and  preventive  repairs  (when  an  inspection 
reveals  the  state  of  the  system  being  worse  than  or  equal  to  k).  The 
convenience  of  a Poisson-type  structure  for  the  number  of  events-per- 
unit-time  made  it  relatively  easy  to  allow  general  freedom  in  the  selec- 
tion of  observation  times. 

The  recursive  equations  of  the  expected  cost  to  the  end  of  the 
cycle,  incurred  subsequent  to  an  inspection  at  whicli  the  system  was 
observed  in  state  i,  are: 


= [Jj  + S Pjjlf  


j=i+l 


: '.■*•11, .. 
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(Vhe  r »■ 


(Xt  )' 

-TFir 


(i  < j < n). 


(4.  2) 


J.  = the  expected  inspection  and  occupancy  costs  to  the 

^ next  event  (an  inspection  or  a malfunction)  incurred 

subsequent  to  an  inspection  at  which  the  system  was 
observed  in  state  i, 

c = the  complementary  probability. 

An  optimum  inspection  procedure  is  a specification  of  the  suc- 
cessive inspection  times  (t/,  0 < i < k}  for  which  the  expected  cost 
is  a minimum.  The  total  expected  cost  per  unit  time  LQ/E(t)  can  be 
optimized  in  terms  of  an  auxiliary  optimization  function  (the  Brender's 
method  [3  ].  [l  ])  •• 


Dq  = Lq  - O'  E(t), 


(4.  3) 


where  E(t)  is  the  expected  cycle  length.  From  (4.3  ) the  expected 
cost  is  transformed  to  the  following  recursive  equations: 


D.  = rj.  -oJ.  + F P..(T.)D.1/Pf-  (t.)  (i=  0,  1 k-1), 

1^1  1 ij  1 j"  11  1 


(4.  4) 


where  J.  is  the  expected  time  elapsing  from  an  inspection  at  which  the 

system  was  observed  in  state  i to  the  next  event  (an  inspection  or  a 

malfunction  detection).  The  minimal  value  of  D^,  for  fixed  k and  a, 

is  obtained  by  finding  recursively  the  i^'s  (i=  0,  1,...,  k-1)  which 
• ❖ 

minimize  D„.  Since  all  the  t.'s  for  j > i are  known  when  D.  is 
0 J 1 

minimized,  t/  is  obtained  by  minimizing  a nonlinear  function  of  a 
single  variable.  By  varying  a and  k,  the  Dq  is  minimized  to  zero 
with  the  optimal  a and  k . The  minimum  total  expected  cost  is  ob- 
tained  and  is  equal  to  O'  , The  optimal  policy  {k  ;t  . ' s(i=  0,  1 , . . . , k - 1 ) } 
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is  als(;  (jbtain«-d. 

The  work  studied  here  is  based  on  a modification  of  the  model 
used  by  Luss.  Our  model  for  deterioration  is  more  general,  but  the 
admissible  strategies  used  here  are  more  restricted.  Here  we  allow 
the  exponentially  distributed  durations  to  have  different  mean  values, 
and  to  be  positively  correlated.  The  joint  distribution  assumed  for 
these  variables  is  introduced  in  Chapter  IL  Some  relationships  of  de- 
pendence among  these  variables  are  studied  in  the  next  section. 

The  presence  of  correlation  between  interval  durations  permits 
the  modeling  of  a rate  of  deterioration  which  can  be  estimated  during 
a particular  realization  of  its  past.  However,  the  lack  of  a Poisson- 
type  of  structure  for  the  events-per-unit-time  makes  it  much  more 
difficult  here  to  allow  general  freedom  in  the  selection  of  observation 
times.  At  present  only  the  simple  case  of  direct  and  instantaneous  ob- 
servation of  deterioration  jumps  has  been  considered. 

Figure  5 shows  a typical  time  history  of  deterioration  and 
replacement.  The  duration  in  state  (i-1),  prior  to  reaching  state  (i), 
is  r^  The  sequence  {r^}  will  be  Markov,  characterized  by  a m.ulti- 

variate  exponential  distribution  described  in  Chapter  II  and  in  the  next 
section.  Reward  functions  will  be  related  to  the  deterioration  state 
and  the  time  spent  in  each  state.  The  decision  rule  specifies  whether 
or  not  to  replace,  when  entering  each  state  i,  on  the  basis  of  the  his- 
tory of  r^  r.  2 The  Markov  property  simplifies  the  decision 

rule  to  a collection  of  C.^  sets  such  that  we  replace  on  entering  state  i 
if  and  only  if  r.  ^ e C.. 

The  objective  is  to  maximize  the  average  reward  per  unit  time: 


▼ 
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= liiii  Tj  (lota)  rewiird  in  (0,  1)) 

i 


K [Duration  between  renewals]  ~ J) 


p.  fl<t:v.ard  per  renewal 


( i.'i  ) 


(•^.6  ) 


(See  Ross  [27]  for  equivalence  of  (^.5  ) and  (4.6  ).  ) The  mean  reward 
per  renewal  is  defined  here  as: 


/e=  E 


N-l  ^i 


S [ c.(t)dt  - 
0 ' 


0 


(4J 


in  which  : 

N = state  at  which  replacement  occurs  (possibly  random). 

Pj^  = replacement  cost  if  replaced  on  entering  state  N 
(possibly  random). 

c.(t)  = reward  rate  when  in  state  i. 

figure  6 shows  several  rev.  ard  rate  time  functions  c(t) 
which  have  been  considered.  When  one  of  these  c(t)  functions  is  spe- 
cified for  a given  problem,  the  c.  (t)  in  ( 4.7  ) are  assigned  values 
f}  .c(t)  with  : 


to  assure  greater  reward  rates  in  less  deteriorated  states. 
The  mean  duration  jf)  in  ( 4.6  ) is  defined  as  : 


(4.8 


E 


L 


N-l 

E r.  + d 

0 ^ 


N 


(4.9 


to  include  a possibly  random  time  d^  for  carrying  out  a replacement 
at  state  N. 


Figure  5: 

History  of  Deterioration  and  Replacement  (n  = 5)  . 


a)  constant 
t 

b)  linear 


c)  constant-after  set  up 


Figure  6: 

Reward  Rate  Time  Functions  . 


While  the  ultimate  objective  is  to  choose  C ^ to  maximize  L,  it 
is  well  known  that  a related  problem  of  maximizing  : 

is  simpler  [ 1 ].  Indeed,  the  which  maximize  L will  be  identical  to 
those  which  maximize  such  that: 

J^^(ot*)  - 0 , where  J!q  (a)  ^ max  jCq  (o)  . (4.1 1 ) 

° {Cj 

Section  4,  3 considers  the  constant  reward  rate  case  in  which  it 
is  found  that  deterioration  rate  information  is  not  useful  (e,  g. , the 
optimal  policy  is  independent  of  the  amount  of  correlation  between  suc- 
cessive state  durationsj. 

Sections  4.4  and  4.5  consider  other  reward  rate  structures  for 
which  the  optimal  policies  do  make  use  of  estimates  of  the  deteriora- 
tion rates  as  well  as  of  observations  of  the  deterioration  level. 

The  next  section  develops  interesting  dependence  properties  of 
the  {r.}  sequence  which  will  be  needed  in  the  optimization  studies. 

4.  2 Relationships  of  Dependence  Among  Multivariate  Exponential 
Variables 

For  easier  and  more  convenient  notations,  we  make  the  following 
changes  in  Equations  (2.  10)  and  (2.  25)  with  means  : 

E[r.]=ni.  (4-12) 


and  correlation  coefficients  : 
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Therefore  we  have 


and 


Vi 


r. 

X 


(i-p)hi_i  (i-p)n. 


2 


nJ  r . , r. 
1-1  1 


f4.14  ) 


f(fO.  rj.  ^2.> 


"n-l)  = 


p 

r 

n-1 

n-2 

(i-p)'"'’  n n 

■ 'o 

i=G  ^ 

i=0 

1- 

■2 \l  T.T. 

i-'i+l  ■ ■ 


•P  Vh;^;  • ■ ' i i"*"! 


exp 


1 


r , n-2  r(l+p) 

0 +^+  Z ^ 


(i-p)\no  >1i.i  i=i 


I;  n->  2. 

(4.15) 

From  our  stationary  correlation  structure,  we  have  a Markov  sequence. 
The  conditional  density  function  is  easily  determined  from  Equation 
( 4.15); 

‘'"ihi.f'i.2 "o'= 

The  Markov  property  will  simplify  replacement  schedule  optimization, 
a procedure  which  will  also  make  use  of  other  properties  of  the  condi- 
tional densities; 


-1 


£(ri|ri-i)=  [ni{l-p)J  exp 


1 


(1-P)\^ni  Hi.i 
For  example,  it  can  be  shown  [23l  : 


‘OU-Pn/ 


(4.  17) 


) 
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E[r.  |ri_jl  = r)j  + '^i-l' 

Var[ri|ri_i]  = r)?[(l-p)^  + 2p(l-p)  ri_j/r?i_l  j . (.4.19) 

These  conditional  moments  show,  e.  g, , that  the  conditional  mean  of 
exceeds  its  mean  in  proportion  to  the  amount  by  which  exceeds  its 

mean. 

Before  deriving  some  dependent  relations  among  r^  and  r^^  we 
review  the  relationships  among  some  notions  of  multivariate  depen- 
dence [ 2 ] . 


Definitions  ( 4.  1 ):  Given  random  variables  S and  T,  we  say  the  fol- 
lowing: (Definitions  5,  4.  1 of  [ 2 ] ) 

(a)  T is  stochastically  increasing  in  S if  ; 

Pf  T > t|S=  s]  , (4.  20) 

is  increasing  in  s for  all  t.  We  write  SI(T/S). 

(b) LetS  and  T have  joint  probability  density  f(s,  t).  Then  f(s,  t) 
is  totally  positive  of  order  Z if  : 


f(Sj,t2) 


fis^.t^) 


> 0 


(4.  21) 


for  all  Sj  < s^,  tj  < t^  in  the  domain  of  S and  T.  We  write  f(s,  t)  is 
TP^  or,  alternately,  TP2(S^  T), 

Theorem  ( 4.  1 ):  TP2(  S,  T)  i=r>  SI(T/S)  (The  orem  5.  4.  2 of  [ 2 ] ) 


Proof:  From  (4,.  21),  we  obtain  : 
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00 

f {(8^,  T )dT 

I 


< 0 


t t 

[f(s^,T)dT  J f(S2,  T)dT 
0 0 


(4.  22) 


for  Sj  < s^.  Adding  the  top  row  to  the  bottom  row  and  converting  to 
ratios,  we  obtain  the  inequality  for  conditional  probabilities; 

P[T  > t|S  = s^]  < P[T  > t(S  = s^] . (4.  23) 

Thus  SI(T/S)  holds.  II 

The  following  theorems  will  help  us  understand  intuitively  the  in- 
formation about  the  rate  of  deterioration  contained  in  the  past  observation. 

Another  statement  in  the  same  spirit  of  Equation  (4.  18)  (i.  e. , that 
a large  r^^  ^ implies  will  tend  to  be  large  ) is  that  r^  is  stochas- 
tically increasing  in  r^  j.  In  order  to  prove  SI(r^  | r.  ^ ),  it  suffices  to 

show  that  f(r.,  r.  .)  satisfies  the  TP.  condition, 

11-1  2 

Theorem  ( 4,  2 );  r^  is  stochastically  increasing  in  r^ 

Proof;  By  substitution  of  Equation  (l-.M)  into  the  determinant  (4,  21) 
and  removing  the  non-negative  factors  it  remains  to  be  shown  that  : 


Det.  = 


'•  0 , 


(4.  24) 


where  (^  = 2 /[  (1  -p)  n/  n£n£_|]  . With  the  power  series  expansion 

of  Bessel  functions,  each  Bessel  function  can  be  expanded  according 
to  the  power  series  : 
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< 4 


z 2k 

(4) 


IU)=  L — , 
“ i._(,  (tir 


Det,  can  be  written  as  a difference  of  products  of  infinite  summations: 


Det.  = 


00 

z: 


: x2i 

(^) 


i=0  (i')‘ 


oc 


2i 


i=0  (il)^ 


2 J 


2j 


oo 

^ 2 
j=0  (jl)^ 


OO 

V* 


2 ' 


j=0  (j!)^ 


00 

= 2 


') 


2i 


i=0  (il)^ 


00 

2 

j=0 


2j 

OO 

- ' S 


C- 


2i 


(jl)^  i=0  j=0  (j!)^ 


For  each  pair  of  non-negative  integers  (i.  j),  f^e  expanded  version  of  Det, 
will  contain  four  summands  with  factors  [il  ji]  (if  i 4 j)>  such 

terms  if  i = j.  If  i = j + m,  m > 0,  then  the  four  terms  can  be  written  as  : 
2(i+j) 

(4.  25) 

The  sign  of  (4.  25)  will  be  determined  by  the  factor  in  brackets 
{ • } . Noting  that  for  0 < s^  < s^,  0 < t^  < t^,  and  m a non-negative 
integer  : 


e 

n^n^a-p)^ 


m ^ m 

<^1  = < <^2=  ^2 


■^1  " S - "^2  " *^2  • 


Since: 


i 


I 


when: 
o-j  < a-^  . 


it  is  clear  that  (4.  25)  is  non-negative.  Finally,  Det.  is  seen  to  be 
non-negative  because  the  infinite  summation  can  be  grouped  into  terms 
like  (4.  25),  or  into  the  corresponding  i = j terms  which  vanish.  This 
completes  the  proof  that  f(r.,  r.  ) is  TP  , which  implies  that  r.  is 
stochastically  increasing  in  r.  j.(| 

The  SI{rJr^  ^)  property  will  be  used  in  the  optimization  deriva- 
tions later,  in  conjunction  with  another  theorem: 

Theorem  ( 4.  3 );  If  SI(x(y)  and  h(x)  is  an  increasing  function  of  x, 
then  E[h(x)  |y]  is  increasing  in  y.  (Proposition  3.  1 on  p.  22  of  f 1 7] . ) 

Proof:  For  Yj  ^ y2  ’ 

ao 

E[h(x)|y^]  - E[h(x)|y2]  = J h(x)  [ f{x  | y^ ) - f (x  | y^)]  dx. 

-00 


But  by  assumption  : 

/[f(x|y,  )-f(x(y2)]  dx  = P[  X > t|yj] 
t 

also: 

00 

/[f(x|yj)  - f(x|y2)]  dx  = 0. 

-00 


- P[x  > tly^]  > 0 
for  all  t. 


These  expressions  allow  us  to  invoke  Lemma  1 and  its  corollary  on 
p,  1 20  of  ( 2 ] , viz.: 

If  W (x)  is  a Lebesque-Stieltjes  measure,  not  necessarily  positive, 


for  which  : 


dW (x)  > 0 for  all  t, 


f dW(x)  = 0. 


and  h(x)  is  increasing,  then  J h(x)dW(x)  > 0.  This  completes  the  proof. 


Another  interesting  property  is  that  the  random  variables 


"O'  "i- 


r^  ^ are  conditionally  increasing  in  sequence.  This  is 


equivalent  to  saying  that  n is  stochastically  increasing  inr^  j,...,  r^,  r^ 


for  i = 1, . . . , n-1,  i.  e. 


Pr  r.  > r.  I r.  , , r.  , 
1 1 ' 1-1  1-: 


T'"ol 


is  increasing  in  r^  f ^i  2’  * * ’ ’ "l’  "o* 


Theorem  ( 4.  4 );  If  the  random  variables  r^,  r^ r^  ^ obey  the 

joint  density  function  in  Equation  ( 4.1 5 ),  then  r^,  r^,...,  r^  ^ are 
stochastically  increasing  in  sequence. 

Proof;  We  invoke  the  Theorem  (5.  4.  14)  of  [2  ],  viz.: 

Let  r^,  ....  r^  ^ have  joint  density  • • • » j) 

which  is  TP^  in  each  pair  of  arguments  for  fixed  values 
of  the  remaining  arguments.  Then  r^,  . . . , r^  ^ are 
conditionally  increasing  in  sequence. 

In  order  to  prove  that  f(r.,  I ~ ^1-2“  ' ' ’ ‘ "o  “ ^0^  ^ ^^2 

^"i’  "i-1^’  suffices  to  show  that  f(r^,  j I 2 “ '^i  2^  ^ "^^2 

(rj.,  because  of  the  Markov  property  of  the  sequence.  We 

have  : 


I 


If 


•39. 


'<V''i.J''i-2=  V2>=  '1-2' 


and : 


= f(r..  r._,) 


"i-2) 


f(Si.ti|ci_2)  f(^.t2|c._2) 


f(s2*  S l‘^i-2^  2’ *2  I *^1-2^ 


-wj 


f(Sj,  tj ) 


f(s2.tj) 


f(sj.t2) 


f(s2.t2) 


which  is  always  nonnegative  from  Theorem  { 4.  2 ),  Therefore  we  have 
^^2  ^^i'  ^1-1^’  those  random  variables  which  are  not  next  to  each 

other,  such  as  r^  and  r^  for  j ^ i-1,  they  are  conditionally  independent 
by  the  Markov  property  of  the  { r^}  sequence.  This  implies  Tp^  (r^,  r^) 
and  completes  the  proof.  1| 

The  stochastically  increasing  property  is  not  sufficient  for  our 
analytical  purposes.  To  show  that  the  conditional  density  function  of 
r^  given  r^_^  is  a convexity-preserving  transformation  [l7],  a property 
which  we  shall  need  in  this  chapter,  one  method  is  to  show  that  TP 

O 

true.  However,  a direct  proof  which  requires  a 3 by  3 ma- 
trix expansion  of  Bessel  functions  is  quite  complicated.  Therefore 
other  methods  are  sought  after.  We  summarize  several  important  facts 
and  theorems  in  the  theory  of  total  positivity  fl  7 1 that  will  be  extensively 
used  later. 
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Definition  ( -4.  2 ):  A real  function  K(x,  y)  of  two  variables  varying  over 
linearly  ordered  sets  X and  Y,  respectively,  is  said  to  be  total ly  posi- 
tive of  order  r (abbr.  TP^)  if  for  all  : 


m 


wc  have  the  inequalities: 


eX,  y.fY,  1<m<r  , 

(4.  26) 


K(x^,yj)  K(Xj,y^)  ...  K(Xj,y^) 

K(x^,  y^)  K(x2,  y2)  •• 


K(x  , y,  ) K(x  , Y 
m ' 1 ' m ^ 2 


K(x  , y ) i 
m 'mi 


(•1.  27) 


If  strict  inequality  holds,  then  we  say  that  K is  strictly  totally  positive 

of  order  r (abbr.  STP  ). 

r 

Remark;  From  the  condition  (4.  26),  we  know  that  if  K(x,  y)  is 

TP  then  it  is  also  O P for  1 < m < r. 
r m — — 

For  illustration,  we  begin  by  citing  one  basic  example[  17].  The  function: 


K(x,  y)  = 


^,0(x)^^(y) 


(4.  28) 


is  TP^(STP^)  on  X X Y,  provided  0(x)  and  1/ (y)  are  increasing  func- 
tions (strictly)  on  the  sets  of  X and  Y,  respectively,  of  the  real  line. 


Lemma  ( 4.  1 ):  If  L(?,  q)  is  TPj.  and  M(q,  C)  is  ^ Ihen  : 

K(P,  O = f n)  M(q,  O da(q)  P e X,  ( e X . 

Y 


is  'iPniin(r,  s)*  (I-emma  3.  1.  1.  (a)  of  fi?].  ) 


1 


-41- 


I 


Theorem  ( 4.5):  If  y)  is  and  (p{x)  and  1/7(7)  are  nonzero  posi- 

tive functions  for  x c X and  y c Y,  respectively,  and  if  L(x,  y)=  K(x,  y) 

(ji  (x)ili{y),  fhfn  I,(x.  y)  in  alno  TP^.  (Theorem  "I,  1.  1.  (a)  of  flVl.  ) 

A direct  specification  of  Lemma  ( 4.  1 ) leads  to  the  conclusion 

that  : 


K(x,y)  = 


/eU(^>"(®)eV(y)P(«)d<r(s), 

-00 


X e X,  y € Y, 


(4.  29) 


is  TP^  provided  u,  V,  o,  and  /3  are  monotonic  increasing  functions,  cr  is  a 
sigma-finite  positive  measure,  and  (4.  29)  exists  absolutely.  (Eq.  (3.  1.9) 
of [17].) 

Theorem  ( 4.  6 );  K(x,  y)  = (xy)  is  TP^  for  0<x,y<oo,  where 
denotes  the  modified  Bessel  ftinction  of  order  a.  (Example  on  p.  1 0 1 of 

[17].) 

Proof;  The  modified  Bessel  function  can  be  expanded  as; 


2n  I 

where  a^=  [2  • ni  . p/o+n+l)]  and  (r(t)  is  a sigma-finite  discrete 

measure  concentrating  mass  at  t = n.  The  integral  has  the  same 
form  as  that  of  Equation  (4.  29).  With  t,  log  x , and  log  y being  mono- 
tonic increasing  functions,  we  have  shown  that  I^  (xy)  is  TP^  from 
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Lemma  ( 4.  1 ) and  Theorem  (4.5). 


Theorem  ( 4.  7 ):  The  joint  density  function  (4.14)  of  r^.  and  r-_  j. 


i = 1.  2 n-l.  is  TP^  for  0 < r.,  r._ 


Proof:  We  have  : 


"i-l 


< » . 


1.  r.)  = 


i-1'  i' ■ (l-p)nj. itli 


r . 

X 


(i-p)ni.i  (i-p)ni 


2 

i-p 


'^’i’^i-i 


''"i-i^ 


and  yC  s/r._,rj  where  ? ? • . is  TP^  from  Theorem  ( 4.  6 ). 


Therefore,  with: 


r. 


T-l 

0-p)»»i_l 


n 


1-1 


and: 


r. 

i 


(1-p  )n: 


Hi 


we  have  : 


f(ri_i.  r.)=  0(r._^)^(r.)lQ(?77-7:)  . 

Moreover,  we  note  that  0(r.  and  i/(r^)  are  nonzero  positive  functions 

for  r^  ^e[0,  oo)  and  r^e[0,  «),  respectively.  Hence  f(rj^_^,  r.)  is  a 

TP  from  Theorem  ( 4.  5 ).  11 
oo 


i 
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i' 


'ihe  conditional  density  function  has  the  same  slructuri- 

as  the  loint  density  function  r,).  It  is  simple  to  prove  th.  iollou- 

inp  theorem: 

Theorem  ( t.  H );  i'he  conditional  density  function  (4.17)  of  r^  on  rj^_j, 
i = 1.  2 n-1  is  TP^  for-  0 r.,  r._,  < .r  . 


Proof:  We  have  : 


f(rjr._j] 


(i-p)n 


p f 


i-] 


(i-p)o.  ■ (i-pTn-7i 


• I 1 r . ~] 
0^  X i-P 


where  C = and  P(f  'r.r.  ,)  is  TP  from  Theorem  ( 4.  6 ). 

l-p  I U ^ 1 -« 


-^i’li-l 


Therefore,  with 


pr 


i-1 


“ (l-p)r)^  \ 1 

^ ''c^''i^=  (TT^ 


(i-p)n, 


we  have ; 

f(rjri_j)  = 0^(r._j)  -'^(r.)  Io(Ur.r._'  ) . 

Since  (b  (r.  ,)  and  I'l  (r.)  are  nonze ro  positive  functions  for  r.  , e[0,  of) 
and  r.  € [0,  rn),  respectively,  we  know  that  f(rj^|r^_j)  is  the  TP  from 
Theorem  ( 4.  5 ).  II 

Remark:  That  f (rjr.  j)  is  I'P^can  also  be  justified  by  simply  observing; 
that  it  comes  from  the  probability  transition  function  of  a diffusion  process. 
See  Karlin  and  McGregor  [ 30,  31 ) . 

The  followirij'  theorem  concerns  convexity-preserving  transforma- 
tions which  we  shall  use  in  analyzing  our  problem. 


♦ 


1 


- ‘I - 


1,3  C < III. st lint  lli-waril  l\att-- 

1 1j.'  . nu.-,!  ,iiii  r.'//<iiil  I.i1<  < use  with  i h pai-1  iinl.i  rly  .siin- 

plf  to  analyzf.  Wo  will  soo  lliaf  as  long  as  I''[  r.  | I J.  ^ 

all  i,  even  if  the  r.  are  not  exponentially  distributed,  the  optimal  rule 
will  be  to  replace  upon  entering  some  critical  state  k^'b  independent  of 


the  observed  durations  r^.  For  simplicity  of  exposition  we  assume  that 
all  p.  = p and  all  d.  = d. 

Based  on  the  problem  statement  in  Section  4.  1 the  optimal  deci- 
sion on  entering  state  j must  maximize  the  mean  lufure  reward  until 
the  next  renewal,  for  a suitable  a.  Here; 


' N-  1 

£.(0')  = E 

B .r.  r.  , 

. . 11  J-  1 

L-J  M 

- aE 

r . r . , 

1 j - 1 


- p-  o d.  (4.  31  ) 


The  optimal  decisions  for  each  state  will  be  found  in  terms  of  O’,  and 
then  the  proper  a*  (for  producing  decisions  which  maximize  L)  is  the 
one  for  which  the  maximum; 


max  Xq(o  ' >=  £ 0^*^ ’ ^ ~ ^ 


(4.  32) 


Optimization  by  dynapiic  programming  begins  with  considering  the 
decision  at  the  last  step,  i,  e.  , on  entering  state  (n-l).  There  are  two 
choices,  to  replace  (R)  or  not  to  replace  (R),  with  corresponding  values  ; 


£ (a;R)  = -p  - o-d. 


(4.  33) 


J!n.i<a';H)=  F(Pn.i''n.ilrn-2^  - o- K [ Cp,!  I J -p-od 

= F;[  (/?p_j-a)rp_j|rp_2]  -p-crd.  (4.  3‘>) 


» 


Clearly,  the  bt  st  <-ie<  ibion  is  not  to  replace  if  and  only  if: 

^ ^ ,(o;R) 

n- 1 ’ n-  d = n- 1 n- 1 


= ^{^n-l'''n-2^  - ° 


(4.  35) 


The  si(jn  of  (4.  3o)  will  be  the  sign  of  non-negativity 

of  all  interval  durations.  Thus  the  best  decision  depends  on  a and 

the  reward  parameters  /j  , but  not  on  the  previously  observed  duration. 

Two  cases  will  be  considered  separately. 

If  (3  a then  the  best  decision  at  stale  (n-l)  is  not  to  replace. 

n-1 

We  will  now  explain  why,  under  this  condition,  it  is  best  not  to  replace 
at  any  slate  less  than  n.  Consider  the  situation  on  entering  (n-2).  We 
have  already  shown  that  it  is  best  not  to  replace  on  entering  (n-l).  Thus 
the  choice  will  be  based  on  a 2 form  : 

Here  we  have  : 


(4.  37) 


by  assumption,  and: 


h.  3«) 


because  all  r^  0 with  probability  one.  Thus  A^_  > 0 for 

all  r >0,  and  it  is  best  not  to  replace  here,  either.  This  argument 
n-  3 

can  be  repeated  for  states  (n-3),  (n-4) 1,  0. 

The  other  case  to  consider  is  < a which  requires  replace- 

ment on  entering  state  (n-l),  if  the  system  ever  reaches  that  slate. 


J 


When  we  cotisidc- r tlie  derision  on  t-nlei'ing  (n-l^),  lln-  is  ■ 


■ 2 n- 


,)= 


n-  l 


■air  .Jr  „ I 
' n- 2 n-i' 


(d-  )'' ) 


which  has  the  sign  of  If  (^^n-2"^^  replacement  is 

optimal  on  entering  (n- 2)  and  (n-2)  is  considered  next.  This  iteration 
may  eventually  reach  a state  (k-l)  where  '*0  it  is  best  not  to 

replace.  Arguments  similar  to  those  for  the  '*  ® case  show  that 

non- replacement  is  optimal  at  all  states  preceding  the  one  which  first 
arises  in  this  backward  iteration  as  a non- replacenient  state. 

In  summary,  in  the  constant  reward  rate  case  £q(C()  is  maximized 
by  a decision  rule  which  says  replace  on  entering  some  state  k < n which 
depends  on  the  reward  parameters  the  a : 


k=  min{i;  (o -(3^)  ■>  0}  . 


Finally,  we  must  choose  o-  so  that  J?  (o*  ) = 0 ,whe re 


k-  I 

- p - od  + O.-ct)  F[r.] 

0 


(4.  40) 


(4.  41) 


Figure  7 shows  a typical  plot  of  <£2(0)  as  a continuous,  piece - 
wise  linear  curve  whose  zero  crossing  (J!q(o'^=)  = 0)  defines  or  and  the 
optimal  replacement  state  k-  for  maximizing  L. 


Example;  Figure  7 shows  that  the  optimal  av'erage  reward  per  unit 

time  is  2-j  when  k-  = 3,  where  = 5,  j = 4,  1^2“  ^3  ~ ^^4  ' *' 

= 0,  p = 5,  d = 1,  n.  = 2 (i  = 0,  I,  2,  3,  4)  and  n = 5,  F rom  Equation 
5 t 

( 4.  41),  we  have  the  following  values  for  Jl!Q(k(t>),  O'),  the  optimal  k is  a 
function  of  a,  which  remains  constant  when  a varies  over  each  inter- 
val < O'  < |3.  : 
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0)  = 25. 


I)  - i-i. 


^q(3.  2)  = 5. 


£q{2.  5)  = -2. 


. 4)  = -7  , 


and  the  optimal  £ ^(k  = = 3,  a ■!<  = 2-^  ) = f> 


4,4  Linear  Revk-ird  J<atc 

The  linear  reward  rate  case  with  c.(t)=2p.t  has  some  very  inler- 
esting  results.  By  utilizing  the  stochastically  increasing  property  of 
• (rf  IV  ^ ),  the  optinial  rule  will  be  seen  to  depend  on  the  observed  dura- 
tions { r^}  , For  simplicity  we  also  assume  that  all  p^  = p and  all  d^  = d. 

Based  on  the  probleni  statement  in  Section  4.  1,  the  optimal  deci- 
sion on  entering  state  j must  maximize  the  mean  future  reward  until 
the  next  renewal,  i,  e.  , £j(o').  For  a suitable  a,  we  have: 


''N-1 

‘ 1 ■'i 


N-1 

r.  , > -aE  I E r. 

li  = J ^ 


^j-l)  -P-Q'd.  (4.42) 


The  optimal  decisions  lor  each  state  will  ])e  found  in  terms  of  Q,  and 
then  the  proper  Q”:'  (for  producing  decisions  which  maximize  L)  is  the 
case  for  which  the  maximum  : 


= -p  - + E 


N-1 

S 


N-1 


3.  r.  - a*  Z r. 
i = 0 ^ ^ i = 0 ^ 


= 0 


(4.  4 3) 


1 


Optimization  by  dynamic  programming  begins  with  considering  the 
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flccision  at  tin-  last  stf-p.  Sinci'  stale  n r«-prebenls  a taiU-fi  (.(u iipoi.t ni, 
we  <!e  1 i iii  t «■  I y replac  *■  it  wlien  it  ent<;rs  state  n.  Next,  we  consider  th» 
decision  to  be  made  on  entering  state  n-1.  There  are  two  choices:  to 
replace  (R)  or  not  to  replace  (H),  with  corresponding  values: 


"^n-l  ' “ P ■ ’ 


- P ' . 


(4.  44) 


(4.  45) 


tor  £.  j (O').  Clearly,  the  best  decision  is  not  to  replace  if  and  only  it 


-^n-l 


= Vl(^n-2-^^n-l-"):i  0 • 


(4.  46) 


Substituting  in  the  conditional  moments  (see  (4.  1 8)  and  (4,  19)).  we  obtain: 


’^n-1 


This  quadratic  function  of  r^  2 can  have  the  possible  con\'<-x  shapes 
shown  in  Figure  8. 

hTxamination  of  (4.47)  and  Figure  8 shows  that  if  A^_j(0)  ' 0, 

then  — A , (0)  ''  0 and  — A , (r  _)  > 0 . This  insures 
dr  n-1  dr  „ n-1  n- 2 

n-  2 n-.2 
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n-  1 ii- 


A 

i 


. r 


^ - r r 

n- A / ^ n-  2 


r 


n-  2 


Figure  8: 

Possible  Convex  Shapes  for  A , (r  .,) 

n- 1 n-  2 


A . crosses  the  r .-axis  ett  most  once.  Therefore 
n- 1 n-2 


3 r ,,  s 0 j A (r  ) > 0 iff  r .,  >•  r 

n-2—  n-t  n-2  n-2  n- 

and  ; 

^n-l<‘^)"f"n-2=V2-'''n-2’  ' 


2 ’ 


(4.  48) 


(4.  4<), 


Considering  the  decision  at  next  earlier  step,  i.  e.  , upon  entering  state 
(n-2),  the  total  reward  until  the  end  of  the  cycle  is  the  sum  of  the  ex- 
pected reward  during  r^  and  the  expected  reward  after  next  transi- 
tion. In  a manner  parallel  to  the  development  of  (4,  46),  we  get: 

R)  - £ .fQ',R)=E{3  .-r^  . - Gr  , J r J 

n-2  n-2  n-2  n-2  n-2  n-2  n-3 


+ I A , (r  .,)f(r  ,,  I r „)d  r 

2.,,  n-1  n-2  ' n-2  n-3  n-2 


n-  2 


= A ( r , 3 Ot) 
n-2  n-3  n-2 


+ I A ,(r  T)U(r  .,-r  .,)f(r  .Jr  .,)d  r . 

Q n-1'  n-2  n-2  n-2  n-2'  n-3  n-2 


= V2<"n-3'^n.2*°')  ■'^n-2(''n-3^  • 


(4.  50) 


-r)2- 


The  first  term  on  tiie  right  hand  side  of  ( 4.  50)  is  written  in  terms  of 

the  similar  quadratic  function  defined  in  (4.  47)  except  for  different 

parameters.  Therefore  it  has  the  same  possible  convex  shapes  shown 

in  Figure  8.  The  second  term  I „(r  „)  on  the  right  hand  side 

n-2  n-3 

of  ( 4.  50)  is  increasing  in  r^  ^ because  of  Theorem  (t.  3) 
and  that : 


A ,(r  U(r  „-r  ' „) 
n-1  n-2  n-2  n-2 


(4.  51) 


is  an  increasing  function  of  r^  Furthermore,  it  is  also  quadratically 


me  reasing. 


Figure  9: 

The  Sum  of  a Quadratic  Function  and  a 
Convex  Increasing  Positive  Function. 


This  resulting  conditional  expectation  I „(r  „)  is  also  convex  in  r 

n-2  n-3  n-3 

because: 


Vl<"n-2)  ^(^n.2-''n.2)  • 
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I 

I 


is  a convex  function  of  r ..  and  Theorem  (4.  10). 

As  a sum  of  a quadratic  convex  function  and  a positive  increasing 

convex  function,  A (r  ) has  the  possible  convex  shapes  as  shown 
n-^  n-3  ^ 

in  Figure  9(c).  The  curves  a and  b of  Figure  9(c)  have  single  crossings  on  r 

It  is  not  clear  whether  curve  c which  has  double  crossings  on  r will 

^ n-3 

1 . # 
happen  or  not.  If  this  does  occur,  we  shall  choose  r „ as  the  cross- 

n-  3 

ing  on  the  right,  therefore  : ‘ 


if 


r r 

n-3  n-3  * 


(4.  52) 


Because  of  the  convexity  of  A „(r  „),  there  are  at  most  two  zero 

^ n-3 

crossings  on  namely,  ^nd  r^  2.  The  global  optimal  region 

for  replacement  is  : 


^^3=  V3  ^ V3  < ^3^  ' 

But  the  interval  {r  : 0 < r < r „}  is  locally  optimal, 
n— o ii— 3 H— 3 

for  |e  I — 0 : 


(4,  53) 


Because 


n-3 


V2'V3>'<V3lV.|<‘‘ V3  . 


r „+  f 
n-3-  ' 


(4.  54) 


Because  of  the  particular  dependence  relation  among  the  r.'s,  the  following 
conjecture  concerning  the  single-aero  crossing  of  A^  ^ seems  to  be  true*: 

*The  conjecture  is  later  found  to  be  correct.  The  proof  and  an  extensive 
study  are  given  in  Appendix  1. 


will 


Conjeclurc:  There  exists  such  that: 


■ '>  i' 


n-.{  n-.i 


and  : 


^"n-:r  ^3  < "n-3^  * 

Such  nice  analytical  properties  have  eluded  us  at  this  time.  How- 
ever, numerical  examples  have  all  demonstrated  that  such  simple  thresh- 
old criteria  for  the  optimal  decisions  are  true.  To  consider  entering 
stales  {n-3),  (n-4),  ....  1,0,  we  have  only  to  repeat  our  argument 
shown  above.  Hence  we  obtain  the  general  recursive  equation  for  the 
decision  function  on  entering  state  (n-i); 


^ • 1.  P .,  O')  = A .(r  . ,,  B ..a) 

n-i  n-i-r  ^n-i’  ' n-i'  n-i-l’  ^n-i’  ' 


^ / Vi+l("n-i'  ^"-i+l’")U(r^_.-r^^_.)f(r^_.|r^_._j)dr^_. 


^n-i("n-i-l'  ^n-i*  + ^n-i^'n-i-l'  ^n-i+l  ' “)* 


(^.55) 


Here  we  have  used  the  same  conjecture'^for  (n-i),  i.  e.: 


2 '’n  i 1 - ® ^ I i)  > 0 ^ ■ I > r ■ , 

n-i-i  n-i  n-i-l  n-i-1  n-i-1  i 


Vi-i  < Vi-i> 


(4.56) 


The  maximal  value  of  £ q{o!)^  for  fixed  ct , is  obtained  by  finding  recur- 
* 

sively  the  r. 's  (i  = 0,  1 n-2)  which  maximize  ^^(a).  Since  all  the 

* * 
r^'s  for  j i are  known  when  A.^j(r., a)  is  to  be  maximized,  r. 

is  obtained  by  searching  a zero  crossing  of  a nonlinear  function  of  a 

single  variable.  Thus  the  model  is  solved  by  the  following  algorithm  : 

2 

See  Appendix  1 , 


rv'Tjti 
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(See  Appendix  I for  the  computational  aspects.  ) 


a.  Guess  an  initial  value  for  or. 

0 * 

b.  Find  ‘£Q(«)t  (i  = 0,1,...,  n- 2)  by  sol  vinf>  recursive 
equations  (4.  55 ). 

c.  If  |j!Q(a')|  < e — 0,  stop  computations  (let  a*  = a),  otherwise, 
estimate  a new  value  for  a and  return  to  b. 

The  following  properties  of  (q)  can  be  used  to  generate  a sequence 
of  cf-values  converging  to  one  which  satisfies  the  inequality  in  c . 

i)  J1!q  (a)  is  monotone  decreasing,  since  j£q  (o)  has  this  property  for  a 

fixed  policy  (See  Eq.  (4.  42));  and  if  (o^) i (oj ) then  the  policy 

used  to  achieve  jCq  (q^)  could  be  used  to  achieve  an  £q(oj)  > (o^  --a  con- 

tradiction , 

ii)  As  shown  in  the  numerical  example  below  (page  58),  a*  (p  = 0)  is  re- 
latively easy  to  compute.  This  value  is  a lower  bound  on  ce*  {p4  0)  . 

Appendix  3 shows  the  flow  chart  of  the  algorithm. 

f-  xamples  (Linear  Reward  Rate  Case):  Tables  3 and  4 present 
results  from  some  numerical  examples  based  on  the  models  in  this  sec- 
tion. Figures  10  and  11  show  reward  and  policies  for  different 
values  of  correlation.  Assuming  that  p = 5.  0,  d = 1 , q.  = 2.  0,  n = 5, 

Pq  - 5>  /^|  = ^2  ~ ^3  ~ ^4  ~ ^ ^5  ~ results  indicate 

that  the  optimal  average  reward  per  unit  time  increases  as  the  correla- 
tion increases.  We  present  details  of  the  computation  involved.  The 
recursive  function  (Equation  (4.55))  is: 

00 

Ajlri.,.  (J;.  »)  = A.(r._|.  (?i,  o)  + ; 0.^,,  »)  f(r.  | rj_|)d  r,  . 


We  list  the  first  term  on  the  right  hand  side  of  the  above  equation  : 
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f 

% 

1 

^.(r._^,l3i.a)=  E{^.rf  - ar.\r._^) 

= + {/3.8p(l-p)  - “P)  r.  + {/3.8(l-p)^  - 2a(\-p)]  , 

(.*.57) 

for  different  i: 

[ U O',  p)  = p^r^  + {8p(l-p)-0'p  } + 8(l-p)^  - 2fl;(l-p  ) , 

! 2,cy,  p)=  2p^r^  +{l6p(Up)-ap}  +16(1-p)2  - 2«{l-p)  . 

I A^(r^,3,a.p)=3p^r'^^  + {24p(l-p)-ap}r^  + 24{l-p)^-2a{\.p)  , 1 

I P)  = ^P^i'o  + {32p(l-p)-fl'p}  + 32(l-p)^  - 2£y(l-p  ) , 

I = 0,  5,  a,  p = 0)  = 40  - 20?  . j 

The  overall  will  be  the  sum  of  and  the  repeated  integral  I.,  The  ' 

optimal  ■£q(o')  for  fixed  value  of  O'  is  : 

oo 

= " p ■ ^ ^0  ’ (4.58) 

’^0 

where  Hr^)  is  an  exponential  density  function  with  the  mean  value  equal 
to  2. 

In  the  following  we  study  the  two  extremal  cases,  i.  e.  , p =1  and 
p = 0.  These  cases  could  be  done  by  analysis  with  minimal  computer 
usage. 

Case  1 (p  = 1);  Ihere  is  only  one  random  variable,  r^,  in  this  case 
with  r^  = ^2~  ^3  " **4  “ ^Q*  optimal  replacement  state  N will  be 

a function  of  r^.  The  equivalent  optimal  average  reward  rate  can  be 
represented  as: 
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I 


li  i^jy;''Q  ~ ^*"0  ■'  0>  '''■’^  replace  on  entering  state  N.  Iherefore,  we  re- 
place on  entering  state  N if  r^  falls  in  the  given  interval; 


N 


a 


a a 

ir,^ 


O'  a 

^2- 


O'  a 


O'  a 


Then,  we  have 


Xo(a)=  -p-ad  (^Qr2.c.rQ)ie  ^ dr^ +/  (0  ^ r ^.ar^)!  e ^ dr^ 
->  , , --2  « , , 


00  ^ - — =. 


(4.60  ) 


a 


_«£££  a a a a 

- -5-0- +40-20' +20' (e  ®+e  ^+e  ^+e  ^)+32e  ®+24e  ^+I6e  '*+8e 


( 


1 


Thf  optimal  Of*  which  satisfies  J!q(Q'*)=  0 is  easy  to  obtain  and  equals 

16  10,  The  optimal  thresholds  {f-  = — . t = 0,  1,  2,  3}  are  obtained 

^ ' i+1 


from  o"t'. 

Case  2 (o  = 0);  There  is  no  dependency  on  the  past.  Therefore  the  opti- 
mal policy  is  to  find  a fixed  strategy  such  as  a critical  state  k*.  We 
definitely  replace  It  on  entering  state  k-^  To  represent  a*  as  a func- 
tion of  p : 


a*(p ) 


N-  1 
E S 
0 


N-  1 

E E r.  + d 
0 ^ 


(4.  61) 


where  N is  a random  variable.  For  the  fixed  k-'-,  we  have; 


a'*(p=0) 


k=:'-l  .. 

E E 3.r. - p 
0 ^ ^ 
k*-l 

E E r.  + d 
0 ^ 


(4,62  ) 


If  we  use  the  same  strategy  for  the  case  when  p ^ 0,  then  it  is  easily 
seen  that  the  average  reward  per  unit  time  is  ^^(p  = 0).  As  the  optimal 
strategy  should  not  work  worse  than  the  above  strategy  which  is  opti- 
mal for  p = 0,  we  have  : 


»*(p)  2.  = ®)  • (4.63) 

With  our  replacement  policy  structure,  the  independent  case  will  have 
the  minimum  average  reward  per  unit  time.  To  find  the  k*  for  the  in- 
pendent  case,  we  have  only  to  look  for  the  maximum  of  the  average 
reward  rate  which  is  a function  of  k. 


» 
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k.- I''ixc-fl  Rcplacem<Tit  Stale Avera);;<-  Ivi-ward  Hairs 

E ( ^ P:rf)  - p 
i = l 
k=  1 

E(  S r ) + d 
i = 0 


5x8-5  , ,2 

2+1  ' 3 


(5+4)8-5  _ 4 

2+2+1 


(5+4+3)8-5  „ 

2+  2+2+1 


Iherelore,  k=I=  = 2 and  we  decide  to  replace  the  part  on  entering  state  2. 


table  3 


Linear  Reward  Rate  Case:  The  Fast  Convergence  of  (or,  £^(0')) 
to  (a*,  J!Q{Qr+')).  Use  the  Linear  Interpolation  Method. 


^ - 5>  Pq  " ^2  " ^3  " ^4  " ^ 

and  Pg  =0. 


First 
Ite  ration 


(13.  40,  0.  0) 

(1  3.  82,  -1. 0981  83) 

; (13.  00,  5,  802) 

I (14.  96,0,  1 15082)  , 
' (16.  10,-0.005)  ' 


Second 

Iteration 


(13.  70,  -0.  528967) 
(15.  00,  -3.  375) 

(1  5.  04,  -0.  237404) 


Third 

Iteration 


(13.  59,  -0.  00181  2) 
(14.  25,  -0.  0688) 

(1  4.  985,  0.  004868) 


table  4 


Linear  Reward  Rate  Case:  The  Optimal  Replacement  Thresholds  . 


p = 5,  d = 1 , n-  = 2,  n : 
3^  = 1,  and  3.  = 0. 

= 5.  3o  = 5,  3 

= 4,  3,  (3 

3=2. 

•J- 

P 

"o 

"2 

0.  0 

0.  0 

4-  00 

+ oo 

+ OO 

0.  25 

0.7 

3.  8 

11.  3 

37.  0 

0.  50 

1.9 

4.  1 

8.7 

22.  5 

0.  75 

3.  0 

4.7 

8.  1 

17.  9 

o 

o 

• 

4.  25 

5.  36 

8.  05 

16.  10 

4.  5 Constant  Reward  Rate  - After  Set-up  (Alignment)  Time 

In  practice,  the  time  of  set-up  after  each  transition  occurs  is 
finite,  and  this  time  should  be  included  in  the  analysis  of  the  model. 

In  fact,  the  inclusion  of  the  set-up  time  generally  causes  major  changes 
in  the  results. 

In  general,  the  set-up  time  will  be  denoted  by  c.  To  begin  with, 
the  reward  rates  after  the  set-up  period  will  each  be  assumed  to  be  con- 
stant. For  simplicity,  we  also  ass\ime  that  all  p^^  = p and  all  d^  = d. 
Similar  to  the  constant  reward  case,  the  decision  rules  for  the  severely 
deteriorated  states  are  to  replace  it  definitely.  Therefore  we  can  find 
an  optimal  critical  state  k*.  By  utilizing  the  stochastically  increasing 
property  of  f(rj|r^  ^ ),  the  optimal  rules  on  entering  state  {i  : i<k'i‘} 
will  depend  on  the  observed  duration  {r.,  0 < j < i-l}  . The  optimal 
threshold  {r^rk*-!  < i <n-2}  will  be  infinite.  We  definitely  replace 


ire  10; 

2ar  Reward  Rate  C 
rching  of  O'*  . 


0 


1=0.  75 


0 


0|  1 2 3 45  67  89  10 

Figure  1 1 : 

Linear  Reward  Rate  Case  , 
Average  Reward  Dependence  on  p . 
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il  on  <•  rit «•  ri iij'  Hiatt-  )<■<“. 

In  this  section  we  will  combine  the  techniques  used  in  Section 
1.  i .xiui  Sect  ion  1 . 1. 

Based  on  the  problem  statement  in  Section  4.  1,  the  optimal  deci- 
sion on  entering  state  j must  maximize  the  mean  future  reward  until 
the  next  renewal  (the  Principal  of  Optimality),  i.  e.  , £ ^(a),  for  a suit- 
able a , we  have  ; 

N-1 

0.(r,-c)u  (r  -c)ir  ] - p - cx  E [ ^ r | r ] - a d . (i.  64) 

i=j  ■'  i=j  J'" 


The  optimal  decisions  on  entering  each  state  will  be  found  in  terms  of  a, 
and  then  the  proper  a-^  (for  producing  decisions  which  maximize  E)  is  the 
one  for  which  the  maximum  : 

i:Q(«*)  = 0 . 


The 


reward  structure  during  each  state  (r.)  for  this  case  has  the  form: 
E[Ci(ri)|ri_j  ] = 0 . P[  r.  < c|r._J  (4.65) 

+ /3.(-c  + E[  rj  r._j,  r.  > c]  ) • P(r . > c | r^_ ^ ) . 


Optimization  by  dynamic  programming  begins  with  considering  the 
decision  at  the  last  step,  i.  e. , when  the  system  enters  state  (n-1). 

There  are  two  choices,  to  replace  (R)  or  not  to  replace  (^),  with  corres- 
ponding values  : 


■^n-l  = - p - Q-d  , 


(4.66  ) 


■ "^•''n-l  l'■n-2^  - P - Q-d  . (4.67  ) 


♦ 


-f,4. 


lot  J'  ,(<•').  tlic  fli  I'lir  rcnc*-: 

/I-  1 


n-l 


v_ 

= I (-  a')r  , f(r  I r „)dr 
■'  ' n-1  n-1'  n- 2 n-1 

0 

00 

+ f (P  i(r  i-c)-ar  ,)f(r  ,|r  „)dr  , 

■J  n-r  n-1  ' n-l  n-1'  n- 2 n-1 

c 

00 

= f (-Q')r  f(r  , I r )dr  , 

J 'n-1  n-1 ' n- 2 n-1 

0 

oo 

+ f P , (r  ,-c)f(r  , I r „)dr  , 

■I  n-1'  M-1  n-1'  n-2  n-1 


0 


oo 

+ f P , r ,f(r  , I r i 

J n-1  n-1  n-1'  n-2  n-1 

c 

X 

- fi  , c r f(r  , I r r,)dr  , 

n-1  J n-1  ' n-  2 n-1 

c 

X 

= f (-Q')r  ,f(r  , I r „)dr  , 

J ' n-1  ' n-1 ' n-  2 n-1 

X 

+ f P , r ,f(r  , I r o)dr  , 

■J  n-1  n-1  ' n-l'  n-2  n-1 

0 

c 

- f (3  r ,f(r  I r „)dr  , 

J n-1  n-1  n-1'  n-2'  n-1 

0 

X 

- (3  , c f f(r  , [ r o)dr  , 

n-1  J ' n-l  ' n- 2 n-1 

c 

= "’n-r") 

0 


% 
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- r /3  r f(r  I r )dr 

J n-1  n-1  n-1'  n-2  n-l 

0 

0-> 

- fi  , c ) i(  r , I r „)dr  , 

n-1  J n-l ' n-2  n-1 


= A (r  ,13  , a) 

n-1  n-2  n-l*  ' 


(4.68  ) 


The  A.,  for  i = 0,  . . . , n-1,  will  have  a form  similar  to  that  for  (n-1) 
except  that  j is  replaced  with  and: 

Ai(ri-i.  ;3i.  a)  = J Gi(ri)f^ri|ri-i)dri  . (4.69) 

0 


The  integrand  of  A.  has  the  poss 
G. 


c 


(aj  /3.  <0 

Figure 

Possible  Shapes  for  1 


e forms  as  shown  in  Figure  1 2. 

G. 


IZ: 

le  Integrand  G^  . 


Going  backto  Eq.  (4. 69),  the  decision  function  will  always  be  negative  if 
j <Qf,  which  requires  replacement  on  entering  state  (n-1)  if  the  sys- 
tem reaches  that  state.  Next  we  consider  the  decision  on  entering  (n-2) 

with  the  assumption  that  S , < o,  then  A „ is  : 

'^n-1  * n-2 
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(4.70  ) 


^n-2<"n-3’^n-2'‘")=  '‘n- 2^^^n- 2 ' "n-3^'^''n-2 

0 

30 

-3n.3/[min(c.rn.2)]f(r„.2kn.3)^"n-2 

“ "^11-3^^11-3’  ^n-2’ 


A (r  ) < 0 for  all  r „ if  ^ < a. 

n-2  n-3  n- 3 '^n-2 


This  iteration  may  eventually  reach  a state  (k-1)  such  that  13,  . > a , 

K — I 

where  (k-1)  e (0,  I,  2,  ...  , n-  2).  We  have  ; 


^k-l^^k-2’  '^k-l’  ""  "^k-l^^k-2’  ^k-l’ 


+ (Future  Reward  Before  Renewal  = 0} 
00 

= (^k-l-">J  ’■k-l^("k-li^k-2^‘^"k-l 
0 

00 

3k.  1 / [min(c.  r^_  ^)]f(r^_  ^ | r^.2)drj^.i 


= / S-l^'-k-l^^^^k-  l^k-2)‘^''k-l 
0 


(4.71  ) 


The  integrand  Gj^  l^*^k  l^  shape  as  shown  in  Figure  12(b). 

which  is  convex  on  rj^  ^ and  increasing  when  1 Furthermore, 

Ak_i  (^k.  2»  ^k- 2' is  convex  on  rj^  2 from  Theorem  (4.10). 

It  is  easy  to  see  that  Aj^  l^^k  2^  increasing  function  for  large 

^k-2'  pose  a conjecture  similar  to  the  one  in  Section  4.  3. 


Conjecture ; There  exists  2 such  that; 


‘^k-l^*'k-2^  ^ ° ^k-2^  ’^k-2  ’ 
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and  th«'  optimal  decision  region  ; 

^'■k-2=  ■■k-2'='-k.2^  • 

We  nov.  consider  the  decision  at  its  previous  state,  i.  e,  , state  (k-2),  the 
total  reward  until  the  end  of  the  cycle  is  the  sum  of  expected  reward  dur- 
ing  rj^  ^ and  expected  reward  after  next  transition.  In  a manner  parallel 
to  the  development  of  (4.70),  we  obtain: 


'^k-2^''k-3’  ^k-2’  ~ ^k-2^^k-3’  \ 2’ 


^k-2’  ^k-1’  ^^^^'■k-2l^k-3^‘^^k-2 

’'k-2 


^k-2^’'k-3^  ^k-2^^k-3^ 


(4.7  2) 


where  Aj^  ^ is  convex  on  all  rj^  ^ and  ir.creasing  over  large  values  of 
The  second  term  (1^^  on  the  right  hand  side  of  Equation  (4.72) 
IS  stochastically  increasing  in  r,  and  convex  on  r,  because  ; 


^k-l<’'k-2^  U(rk-2-^k-2^ 


is  an  increasing  and  convex  function  of  rj^  2 from  Theorem  ( 4.  3 ) and 
Theorem  (4.  1 0). 

As  a sum  of  a convex  function  and  a positive  increasing  convex 

function,  Aj^_2(rj^  is  also  a convex  function  of  rj^  g.  We  use  the 
3 

same  conjecture  : 


"k.3  ^ ° ^ ^k-2<''k-3)  ^ 


0 iff 


’•k-3>  ’■k-3 


and : 


See  Appendix  1. 


L 
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I 

I 


1 


f 


k-.r  k-3 


\ 

( 


When  considoring  entering  states  {k-3),  (k-*!) (k-i),  ,..,1, 

e repeat  our  argument  shov.n  above.  Wo  forni  the  general  recursive 
equation  for  the  decision  fund. on  on  entering  state  (k-i)  ; 


■^k-i^^k-i-r  '\-i’  ■ \-i^'‘k-i-l’  ^k-i’ 


^ 1 \-i+I<"k-i-\-i+l'^>f("k-ii"k-i-l)^"k-i' 


k-i 


(4.73) 


Here  we  use  the  same  conjecture  for  k-i,  i.  e.  : 


a 


r . , > 0 
k-i-1 


'^k-i^^'k-i-l^ 


0 if 


^k-i-1  ^ ■'k-i-l 


and  : 


C,  Aa)  - [r  , , : 

k-i  k-i-1 


r,  . , < r, 
k-i-1  k 


-.-1^ 


The  maximal  value  of  £^{0),  for  fixed  a,  is  obtained  by  finding  recur- 
sivcly  the  r^  (i  = 0,  1,  , k,  . . . , n-2)  which  maximize  £^(a).  Since 

i*,' 

all  the  r.'s  for  j > i are  known  when  A.,,  is  to  be  maximized,  is 

J 1+1  1 

obtained  by  searching  a zero  crossing  of  a nonlinear  function  of  a single 
variable.  Thus,  the  model  is  solved  by  the  following  algorithm:  (see 
Appendix  4 for  the  flow  chart) 

a.  Guess  an  initial  value  for  O'  (possibly  O'--'  for  constant 
reward  case  with  the  same  reward  structure). 

b.  Find  k,  such  that  /3.  < O'  (i  = k,  . . . , n-1).  Let 
r,  = 00  (i  = k-1,  . . . , n-2). 

c.  Find  ■^q(o')  = max{ £ } , r . (i  = 0,  1,  ... , k-1)  by  solving 
recursive  equations  (4.73). 
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Q 

cl.  If  !l'Q((/)|<r  *0,  flop  computations,  (let  k*  = k,  O"- = <y) 

othcruise,  estimate  a nev,  value  for  a and  return  to  b, 

■^.6  Conclusions 

The  optimal  replacement  of  deteriorating  parts  based  on  the  de- 
terioration process  described  by  a multivariate  exponential  distribution 
was  examined.  For  the  constant  reward  rate  case,  the  optimization  pro- 
cedure does  not  benefit  from  the  history  of  deterioration.  However,  we 
do  utilize  this  knowledge  to  optimize  the  linear  reward  rate  case  and  the 
constant  reward  rate  after  set-up  time  case. 

Nunierical  examples  show  that  the  average  reward  per  unit  time 
for  certain  reward  structures  is  an  increasing  function  of  the  correla- 
tion p.  Due  to  our  assumption  that  the  transition  times  are  directly 
monitored,  the  prediction  of  the  future  after  each  transition  should  de- 
pend on  the  history  and  the  correlation  between  the  past  and  the  future. 
The  theory  we  suggest  is  that  if  the  correlation  increases,  then  the 
certainty  of  the  future  also  increases.  Therefore,  the  optimal  policy 
is  more  efficient  as  the  certainty  of  the  future  increases. 
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roNCLusroNS 

M nil  I va  r lali'  cxpoiifnl  lal  (hslribulmns  havi-  bcj’ii  iiit.ioduci:<l  into 
both  downtime  modeling  and  deterioration  modeling.  The  one  based  on 
the  sum  of  the  squares  of  multinormal  variables  is  selected  here, 
because  of  its  physical  motivation  and  analytical  simplicity. 

The  sums  of  multivariate  exponential  variables  are  found  to  well 
approximate  the  frequently  used  lognormal  downtime  distribution,  ex- 

' 

cept  for  the  third  and  fourth  moments.  If  the  subsidiary  times  of  a down- 
time can  be  characterized  as  exponential  variables,  then  the  lognormal 
distribution,  to  our  Judgment,  is  used  just  for  the  purpose  of  expediency. 

Several  dependence  properties  are  developed  for  this  particular 
kind  of  distribution.  The  monotonicity-preserving  transformation  was 
applied  by  Barlow  and  Proschanf2]  to  the  reliability  theory.  Here  we 
have  introduced  the  convexity-preserving  transformation  to  problems 
in  the  same  area. 

Different  reward  structures  such  as  constant,  linear,  and  set-up 
cases  are  considered  for  deterioration  modeling.  Only  incomplete  re- 
sults have  been  found  concerning  the  conjecture  that  there  is  at  most 
one  zero  crossing  of  the  decision  function.  However,  we  have  shown 
that  there  are  at  most  two  zero  crossings.  Also  the  numerical  exam- 
ples have  all  demonstrated  that  we  have  only  single  decision  threshold. 

Areas  for  further  research  include  both  extensions  of  our  results 

I 

to  gamma  marginal  distributions  and  nonstationary  correlations.  j 

It  could  also  be  interesting  to  consider  systems  of  more  than  one 
component,  i.  e.  , correlations  exist  both  in  parallel  and  in  series. 
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I 


The  sclu*f]ulin(4  of  inspections  for  systems  in  v.hirh  de^'ree  of  de- 
terioration can  be  observed  only  through  inspections  is  difficult  but 
t hallenging. 

We  modeled  downtime  and  deterioration  separately.  It  is  interest- 
ing to  consider  them  together.  The  time  to  carry  out  replacement  of  the 
deterioration  modeling  can  be  a downtime  process. 


1 

j 
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APPENDIX  1 


PROOF  OF  THE  ZERO-CROSSING  THEOREM 


The  conjecture  concerning  the  zero-crossing  properties  of  the  decision 
functions  in  Chapter  4 is  proved  here  in  a more  general  form.  Let  the 
durations  in  states  0,  n-1  be  denoted  by  r^,  r^ r^  ^ and  as- 
sume that  r.,  0,  5.15.*’'^'  ^ Markovian  sequence  with  stationary 

probability  transition  function: 

f_  (yix)  = f (y|x)  = f (y|x) 

^ i-1  10 


for  all  2 < i < n-  1,  where  f (y  I x)  is  given  in  Equation  (4.  17);  note  that 

'‘l  I *^0 

Hj  is  not  necessarily  equal  to  As  we  have  already  shown  this  conditional 

density  function  is  TP^  and  it  preserves  convexity  when  used  as  a kernel  in 
an  integral  transformation.  We  also  assume  that  at  state  i the  replacement 
cost  and  time,  i.  e.  , p^  and  d^  are  non  negative  random  variables  independent  of 
r^,  0 J i-  1.  The  reward  rate  at  state  i is  again  assumed  to  be  equal  to 
c (t)  where  c (t)  is  a nonnegative  function  and  jS^  ^ ^^0  for  0 i ^ n-  1 , 

= 0.  The  problem  is  to  find  the  optimal  stopping  time  N which,  for  a 
fixed  O'  > 0,  .Tiaximizes: 


N-1  “"i  N-I 

E [ L ^ - Pn  ■ " 1 > 

i=0  ^0  IN  IN  ._Q  i 


(Al.  1) 


whe  re  N depends  only  onr^,  0<^i^N-l.  Because  of  the  Markov  property, 
an  equivalent  problem  is  to  find  the  optimal  decision  rule  which  requires 
that  the  part  be  replaced  upon  entering  state  j iff  r^  ^ € C.  (a).  In  this 
Appendix  we  give  a sufficient  condition  under  which  c.{q)  has  a very  simple 
form,  namely  ^^(0)  = [ 0,  r ^ j ).  We  need  the  following  definitions  and 


?» 
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theorem  about  the  variation-diminishing  property  of  total  positive  kernels. 

All  the  material  is  taken  from  Karlin  [ 177  with  only  slight  changes. 

Let  h(t)  be  a function  defined  in  [ 0,»)  and  let: 

S'(h}  = S'fh(t)J  = sup  S'  [h(tj),  hlt^) ^ t 

arbitrary  but  finite,  and  S'  (x^,  x^,  ....  x^)  is  the  numoer  of  sign  changes 
of  the  indicated  sequence,  zero  terms  being  discarded. 

Let  k(x,y)  defined  on  [0,oo)  x [0,«5)  be  Borel-measurable  and  assume 

30 

that  the  integral  Jk(x,  y)dp(y)  exists  for  any  x in  [ 0,  »).  Here  M-  re- 
0 

presents  a fixed  sigma-finite  regular  measure  defined  on  [ 0,  »)  such  that 
p.(U)>  0 for  each  open  set  U for  which  U'^[0,oo)  is  nonempty.  Let  h be  bounded 
and  Borel-measurable  on  [0,'X)),  and  consider  the  transformation: 

00 

g(x)  = (Th)(x)  = / k (x,y)  h (y)d  p (y) 

0 

Theorem  (A  1.  1)  (Theorem  3.  1 on  p.  21  of  Karlin  [ 17]): 

If  k is  TP^  and  satisfies  the  integrability  requirements  stated  above,  then: 

S'(g)  = s'  (T  f)  < S (h)  provided  S (h)  < r-  1 . 

For  the  case  in  which  h is  piecewlse-continuous,  if  S (g)  = S (h)  < r-1, 
we  further  assert  that  the  values  of  the  functions  h and  g exhibit  the  same 
sequence  of  signs  when  their  respective  arguments  traverse  the  half  line 
[ 0,  oo)  from  left  to  right. 

Theorem  (A  1.  2)  Let  h(y)  be  a Borel-measurable  function  such  that 
|h(y)|  < a + by  for  some  a > 0,  b > 0,  and  positive  integer  m.  Let: 

CO 

g(x)  = / h(y)  f (y|x)  dy 

0 

where  f (ylx)  is  the  conditional  density  function  of  r^  given 


Then: 


-74- 


(i)  S‘(^)  V S'  (h), 

(ii)  there  exists  a'  and  b'  such  that  | g(x)  a'  + b'  x 

Proof:  We  can  write: 

oci 

g(x)  = f [h(y)  / (a  + b y f (y|x)  d^(y) 

0 

where: 

d^  (y)  = (a+  by*^*^  d y. 

It  is  clear  that: 

00 

f f (ylx)  dp.  (y)  , 

0 

^ 2m 

exists  for  all  x e [ 0,  ■«)  and  that  h(y)  / (a  + b y ) is  bounded.  Since  f (y  ] x) 

is  TP^  , by  the  above  theorem  S (g)  < S (h).  Next  consider: 

1 I g(x)  1<  / |h(y)  I f(y  |x)  d X < a + b J y^”^  f(y  |x)  dx  , 

* 0 0 

00  2m 

but  it  is  a well  known  fact  that  f Y ^ f(ylx)dy  is  a polynomial  of 

0 

degree  at  most  2m  [33],  therefore  (ii)  is  established. 

Going  back  to  our  problem  oi  (Al.  1).  If  we  assume  the  reward  rate 
function  c(t)  is  a non-decreasing  function  not  growing  faster  than  a certain 
polynomial  (e.  g.  , a constant,  a linear  function  in  time  or  either  of  the  two 
! after  set-up).  Then: 

K 

w(t)^  J c(t)  dt  , 

i ^ 

is  a continuous  convex  function,  w (0)  = 0,  and  growing  not  faster  than  a 
certain  polynomial.  We  shall  also  write  E(p^+  od^^)  = e^(o)  for  1 <_  i £ n. 

By  using  the  dynamic  programming  approach  to  solving  this  optimiza- 
tion problem  we  have  the  following: 


V 
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•^n-  1 1 ■ " ^n-  J ^n-z'"  "^n-  1 ’ 

^n-I  -‘^^n-1  ' '"n  I ^^n-Z  1 

=^n-l("'"n-2^  = ^ [ ^^n- 1 1 > "n- 1 ' "n- 2 ^ * 

and  in  general  , 1 ^ j < n-Z: 
jCj  (a,  R)  = -e^(Q')  , 

f.(a,R)  = E [^jW(r.)  - or. 

+ max  (X.^jla.R),  (a,  R))]  rj_j  ] 

= E[3.w(r^)  - or.  - e._^^(o}  +(^j+i  (o’))'*’  Uj.!  ] . 

where  for  1 <.  j <^  n-  1; 


A .(o)  = A . (o,  r^_j)  = £.  (o,R)  - £.  (a.R) 

= E[0.w(r.)  - ar.+  (A  (a))'*’ I rj_  j ] + e.(a)  - e._^^(Q')  . 

and  (x)^  = max  (0,x)  . Then  the  optimal  decision  rule  upon  entering 
state  j is:  (i)  to  replace  the  part  if  A^(q',  j)  < 0;  (ii)  not  to  replace  if 
A j (o,  r.  j)  > 0;  (iii)  either  one  if  A^  (a,  r^  j)  = 0 . 


Theorem  (A  1.  3):  Let  o be  fixed,  o >0  . If  the  following  assumption  are 
satisfied: 


(i) 

(ii) 

(iii) 


• (*„=»• 

c(t)  is  non-decreasing  and  bounded  above  by  a polynomial, 
the  sequence  e.(a),  is  convex  and  non-decreasing. 


i.  e.  , 
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e < e.  , and  e.  < (e.  ,+  e.^,)/2,  then  there  exist  0<r  Ir,!--- 

j - j+l  J - J- 1 J+^  u 1 

< < oc  such  that  the  optimal  decision  rule  upon  entering  state  j is: 

— n-  1 — 

3{C  ^ 

(i)  to  replace  if  r^  I ^ 1 ‘ not  to  replace  if  rj_^  — ^j-1  ' 

Proof;  Starting  with  the  (n-1)  th  state,  the  difference  in  expected  reward 
of  replacing  and  not  replacing  the  part  given  past  observation  is: 

'^n-l  ^ ^n-1  ^‘^’''n-l^l^n-2^  ’ 

where: 

^n-l<  “'^n-l>  « ‘‘n- 1 ^n-1  ' 

Since  w(r  , ) is  a continuous  convex  function,  thus  G , (a.  r ,)  is  also 
n- 1 n- 1 n- 1 

continuous,  convex,  bounded  in  absolute  value  by  a polynomial,  and  by 
assumption  (iii): 

^n-  1 ~ ®n-  1 ' 

This  implies  that  G , (p,  r)  can  have  at  most  one  sign  change,  i.  e.  : 

^ n- 1 

S"  (G  , (a, r ,))<!. 
n- 1 n- 1 — 

We  will  now  prove  the  following  lemma  which  we  shall  need  later. 

Liemma  (A  1.  1):  If  h(y)  is  ccntinuous,  convex,  bounded  in  absolute  value  by 
a polynomial,  and  h(0)  ^ 0,  then: 


30 

I 


h(y)  f (y|x)  dy  , 


is  also  continuous,  convex,  bounded  in  absolute  value  by  a polynomial,  and 
must  belong  to  one  of  the  following  three  categories: 

(I)  g(x)  > 0 for  all  X > 0 , 

(II)  g(x)  < 0 for  all  X > 0 except  with  a possible  zero  at  x = 0, 


-77- 


3|t  Jjt 

(III)  there  exists  a unique  x , 0 < x < ^ , such  that  g(x)  > 0 

for  all  X > X and  g (x)  < 0 for  x < x except  for  a possible 
zero  at  X = 0 . 

Proof;  That  g(x)  is  continuous  can  be  verified  by  examining  f(ylx)  of 
Eq.  (4.  17)  directly.  Since  h(0)  < 0,  v,e  have  S (h)  < 1,  thus  by  Theorems 
(4.  9)  and  (Al.  2)  g(x)  is  convex,  bounded  in  absolute  value  by  a polynomial, 

_ - 5jC 

and  S (g(x))  £ 1.  If  S (g)  = 1,  then  there  exists  a unique  x where  the  sign 

changes,  furthermore,  S (h)  must  be  equal  to  1.  h(0)  < 0 implies  that 

the  sign  of  h(x)  is  first  negative  then  becomes  positive  as  x goes  from  0 to 

oo  . By  the  second  half  of  Theorem  (Al.  1)  this  holds  for  g(x)  also,  thus 

g(x)  is  in  the  third  category.  If  S (g)  = 0,  then  either  g(x)  > 0 for  all  x > 0 

which  is  in  category  1 or  g(x)  ^ 0 for  all  x ^ 0 and  with  strict  inequality 

for  some  x;  by  convexity,  g(x)  is  in  category  II.  || 

We  will  now  resume  the  proof  of  the  theorem.  Since  G , (a,  r) 

n-  i 

satisfies  all  the  requirements  in  the  above  lemma,  thus  A , (a,  r _) 

^ n-1  n-2 

belongs  to  one  of  the  three  categories  defined  in  the  lemma.  Therefore  we 

>1=  sf:  sic 

can  choose  for  each  category,  r _ where:  (I)  r _ = 0:  (II)  r ^ = oo  : 

^ ' n-2  n-2  n-2 

5|C  S|C 

and  (III)  ^ = -K  . It  is  clear  that  this  decision  rule  is  optimal  because 

for  all  r > 0: 
n-2  — 


1 ^ ® (or  < 0 ) 

n-1  n-2  — — 


if.  r > r*  (or  < r*  ) . 
n-2  — n-2  n-2 


At  state  n-2,  we  have: 


^{a)  = , w (r 

n-2  •^i^n-2  n- 


2'  - O' 


^<"n-2l"n-3)‘^"n-2 


and: 


I 
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- ^n-1 

00 

= - ^n-l>  + «n-2  ^ e^_  ^ {«) 

+ e^(Q')  ] f (t  j r ) dt  . 

With  the  assumption  that  ^ ^ + e__(o)  > 0 , 

thenA^_^  ia,  r)  - (a,  r)  >0  for  all  r. 

Similarly,  for  l<jln-3: 

- ^j+i(a')  = ) w (t)  + (A.^^(o))+ 

- +ej  (O')  - 2e._^j(a)  + ^(t|r)  dt  , 

therefore 

A^  (a,  r)  > (a,  r)  for  all  r > 0 if 

^j+1  («.«•)>  all  r > 0 . By  mathematical  induction  we 

have  proved  the  following  lemma, 

L_emma  (Al.2)^  Under  the  assumptions  of  Theorem  (Al.  3) 

A^  (o,  r)  > A^_^^  (a,  r).  for  ail  r > 0 and  for  1 < j < n-2. 

The  relation  between  A . (a)  and  A.^j(a)  is  summarized  in  the  following 
lemma: 

I^mma  (A  I.  3);  If  A.^^  (q,.  r .)  possesses  all  the  properties  of  g(x)  in 

Lemma  (Al.  1),  then  so  does  A (a.  r ) 

J ’ j-r- 

Proof;  By  definition: 

00 

lo'.-l  = ^OjW(t)  -<rt  + (A,^,(o.t))%a.w  -e.^,  to)]  £(t|r|  dt. 

(Al.  2) 


V 
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SLnce  ^ assumed  to  be  convex,  so  are  (A^^^(a,t))  and  the 

integrand  of  Eq.  (Al.  2),  thus  (a  , r)  is  convex.  In  addition,  A^(a,  r)  is 
continuous,  bounded  in  absolute  value  by  a polynomial. 

If  A^^^(q')  is  in  category  I,  i.  e.  , A^^^  (o-)  > 0 , then  by  Lemma  (A  1.  2) 

A.  (a)  > A.^^  (a)  > 0 is  also  in  this  category.  If  in  category  11 

or  m,  then  A.^j(a,  0)<0.  Therefore  the  integrand  of  Eq.  (A  1 . 2)  at  t = 0 

is  equal  to: 

e.{a)  - e.  , (a)  < 0 , 

J J+1 

by  assumption  (ii),  thus  the  integrand  satisfies  all  the  requirements  on 

h{y)  in  Lemma  (Al.  1).  By  the  same  lemma,  A.  (a,  r^_^)  has  all  the  properties 

of  g(x).  II 

Since  A , (a)  possesses  all  the  properties  of  g(x)  in  Lemma  (Al.  1), 

A.(o)  also  has  all  these  properties  by  the  above  lemma.  We  choose  r^ 

^ ^ — n 

according  to  the  following  rule:  If  A^  (a)  is  in  category  I then  r^_  ^ 

If  A.  (o)  is  in  category  II  then  r ^ = <»  • If  Aj  (a)  is  in  category  III  then 

r*  is  chosen  as  the  unique  zero-crossing  point  where  A (o,  r)>  0 for 
J-1 

r > r*  J . It  is  obvious  that  the  decision  rule  based  on  I'j.j  optimal  at 
state  j;  thus  r*  for  0 < j < n-2,  gives  us  a decision  rule  which  is  optimal. 

^ ^ iji 

By  Lemma  (Al.  2),  if  r . = 0 (i.  e.  , category  I)  then  r"  j = 0 ; if  0 < r . < * 

J 

(i.  e.  , category  111)  then  by  the  fact  that  A (a,  r)  > A^^^(q',  r)  > 0 for 

all  r > r*we  have  r*  , < r*  . If  r*  = <»  , then  obviously  r ^ ^ < r . This 

j J-^JJ  jj 

completes  the  proof  of  the  theorem.  || 


Corollary:  If  E(p.)  = E(p^)  and  E (d^)  = E(d^)  for  1<  i < n , then 

eAa)  = e (a)  for  1 < i < n and  the  theorem  holds, 
i n ” 

Remark  1:  Theorem  (Al.  3)  holds  for  the  special  cases  discussed  in 
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CyliapltT  I with  conm.tiU  rrijlHcriiioiU  cunt  liiidt  time  111 nntgliuul  all  llur  cJil- 
fr>rr>n1  ofAtf*. 

Remark  Z:  In  the  proof  of  Theorem  (Al.  3)  we  "have  relied  heavily  on  the 
stationarity  of  the  probability  transition  function  and  assumption  (i)  and  (iii) 
to  keep  the  number  of  sign  changes  of  the  decision  functions  not  greater  than 
one.  However,  if  the  conditional  density  is  allowed  to  be  nonstationary  and 
the  only  requirement  is  that  c(t)  is  non-decreasing,  discarding  assumptions 
(i)  and  (iii), then  (q,)  is  still  convex  but  may  have  two  nontrivial  zero- 

crossing points. 


-81- 


APPENDIX  2 


COMPUTATIONAL  ASPECTS  OF  LINEAR 
REWARD  RATE  CASE 


In  our  computation  of  integrals,  the  interval  for  "Simpson's  Rule" 
is  0.  1.  The  lowest  sUndard  deviation  we  have  considered  is  0.  5.  There- 
fore,we  have  a good  approximation  to  the  integral  by  the  "Simpson's  Rule." 

The  memory  requirements  are  low.  The  minimum  86k  memory 
is  enough.  But  the  integrations  are  very  time  consuming,  because  there 
are  a few  nested  DO-loops  that  we  have  to  carry  out.  For  each  cycle  in  a 
nested  DO-loop,  the  main  program  has  to  call  the  Bessel  function  sub- 
routine. This  subroutine  takes  about  one  millisecond  before  it  returns 
to  the  main  program.  Therefore,  the  computation  is  very  costly  even 
for  a fixed  value  of  a . We  cannot  put  another  DO-loop  for  a . 

The  method  of  search  for  is  to  guess  an  initial  value  for  a-  (say  a^)  , 
then  the  second  value  for  a*  (say  "'ill  depend  on  the  sign  of 
Hence  we  have  ■A)  (oj)  and  we  can  use  linear  interpolation  to  get 

a good  approximation  of  a*.  The  following  is  a FORTRAN  program  for  our 
problem,  where  n = 5,  = 5,  = 4,  p^  =3,  p^  = 2,  p^  = 1,  p^  = 0 , 

C .(+)  = 2p.t  . 0 < p < 1 . 

DIMENSION  D4(200),  D3(200),  D2(200),  Dl(200) 

DIMENSION  Z4(200).  Z3,  (200),  Z2(200),  Zl(200) 

DIMENSION  V(200),  Y(200) 

D = l.  0 

C THE  REPLACEMENT  DURATION 

P = 5.  0 

C THE  REPLACEMENT  COST 

R = 0.  NNNN 

C THE  CORRELATION  ROO 

F = NN.  NNNNNN 

C THE  BRENDER'S  PARAMETER  "ALPHA" 

A = 1.  /(2.  -2.*R) 

B = A 
C = R#A 

E = SQRT(R)/(1.  -R) 
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C THE  PARAMETERS  FOR  THE  CONDITIONAL  DISTRIBUTION 
G = 2.  0*F*(1. 0-R) 

H=8.  0*(l.  0-R)*-2 
Q=8.  0*R*(1. 0-R) 

T = R*F 
S=R-*2 

C THE  COEFFICIENTS  OF  REWARD  DECISION  FUNCTIONS 

R3S=(-(Q-T)+SQRT((Q-T)-v2-4.  0*S*(H-G)))/(2.  0*S) 

KK3  = (10,  0)-R3S-l.  0 
C FIRST  ZERO  CROSSING 

DO  15  1=1. 200 
V(I)=0.  l-(I-l) 

DO  5 K=l,  200 

D4(K)=0.  01*S=:'(K+KK3)--2+(Q-T)*(K+KK3)*0.  1+H-G 
W4=A*((K+KK3)*0.  1) 

X4=E*SQRT(0.  01*(I-1)*(K+KK3)) 

CALL  BESI(X4,  0.  BI,  lER) 

Y(K)=D4(K)*A*EXP(-W4)*BI 
5 CONTINUE 
SIMP=0.  0 
DO  10  M = l,  197.  2 

SIMP=SIMP+(0.  1 /3)*(Y(M)+4.*Y(M+l)+Y(M+2)) 

10  CONTINUE 

Z 3(I)=SIMP*EXP{-((I-1  )*C*0.  1)) 

D3(I)=0.  0 2*S*(I-l)**2+(2.  0*Q-T)*(I-1)*0.  1+2.  0*H-G+Z3(I) 

15  CONTINUE 
DO  16  1=1, 200 

IF  (D3(I))  16.17.17 
17  KK2=I 

C SECOND  ZERO  CROSSING 
1=200 

16  CONTINUE 
DO  35  1=1,  200 
DO  25  K=KK2.  200 
W3=A*K*0.  1 

X3=E*SQRT(0.  01*(I-1)=:=K) 

CALL  BESI(X3,  O.BI.IER) 

Y(K)=D3(K)*EXP(-W3)*BI 
25  CONTINUE 
SIMP=0.  0 

DO  24  M=KK2,  197.  2 

SIMP=SIMP+(0.  l/3)*(Y(M)+4.  *Y(M+1  )+Y(M+2)) 

24  CONTINUE 

Z2(I)=SIMP*EXP(-((I-1  )*C*9.  1))*A 

D2(I)=0.  0 3*S-(I-l)**2+(3.  0*Q-T)*(I-1)*0.  1+3.  0*H-G+Z2(I) 

35  CONTINUE 
DO  36  1 = 1,  200 

IF  (D2(I))  36.  37.  37 
37  KK1=I 
1=200 

36  CONTINUE 


. < 
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DO  55  1=1, 200 
DO  45  K=KKt . 200 
W2=A*K*0.  1 

X2=E*SQRT(0.  01*(I-1)*K) 

CALL  BESI(X2.0.BI,IER) 

Y{K)=D2{K)*EXP(-W2)*BI 
45  CONTINUE 
SIMP=0.  0 

DO  44  M=KK1.  197.  2 

SIMP=SIMP+(0.  l/3)*(Y(M)+4.  * Y (M+1 ) +Y (M+2)) 

44  CONTINUE 

Z1  (I)=SIMP-EXP(-((I-1)-C*0.  1))=^A 

D1(I)=0.  04*S*(I-l)**2+(4.  0*Q-T)*(I-1)*0.  1+4.  0*H-G+Z1(I) 

55  CONTINUE 
DO  56  1=1.  200 

IF  (Dl(I))  56.  57.  57 
57  KKZ=I 
1=200 

56  CONTINUE 

DO  65  K=KKZ.  200 
W1=0.  05+K 

Y{K)=D1(K)*EXP(-W1) 

65  CONTINUE 
SIMP=0.  0 

DO  64  M=KKZ,  197.  2 

SIMP=SIMP+(0.  l/3)*(Y(M)+4.  *Y(M+1  )+Y (M+2)) 

64  CONTINUE 
ZZ=SIMP-0.  5 
RW=-P-F*D+40-2,  0+F+ZZ 
WRITE  (6,  101)(RW,KK  + .KK2.KK1,KKZ) 

101  FORMAT  (5X,  FIO.  6,  5X.I4,  5X.I4.  5X,I4.  5X.I4) 

WRITE(6 . 1 1 1 )(V(I).  Z 3(1) , D3(I) . Z 2(1) . D2(I) . Z 1 (I) . D1  (I) , 1=1 , 200) 

111  FORMAT(5X.  F4.  1,  3X.  FIO.  4,  3X,  F9.  4.  3X.  F9.  4.  3X.  F9  4.  3X.  F9.  4.  3X.  F9.  4) 
STOP 
END 
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I 

APPENDIX  4 

FLOW  CHART  OF  CONSTANT  REWARD  RATE 
AFTER  SET-UP  TIME  CASE 


k 
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